{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ce02ad0c",
   "metadata": {},
   "source": [
    "# Drift detection on molecular graphs\n",
    "\n",
    "## Methods\n",
    "\n",
    "We illustrate drift detection on molecular graphs using a variety of detectors:\n",
    "\n",
    "- [Kolmogorov-Smirnov detector](https://docs.seldon.io/projects/alibi-detect/en/stable/cd/methods/ksdrift.html) on the output of the binary classification [Graph Isomorphism Network](https://arxiv.org/abs/1810.00826) to detect **prediction distribution shift**.\n",
    "\n",
    "- [Model Uncertainty detector](https://docs.seldon.io/projects/alibi-detect/en/stable/cd/methods/modeluncdrift.html) which leverages a measure of uncertainty on the model predictions (in this case [MC dropout](https://arxiv.org/abs/1506.02142)) to detect drift which could lead to **degradation of model performance**.\n",
    "\n",
    "- [Maximum Mean Discrepancy detector](https://docs.seldon.io/projects/alibi-detect/en/stable/cd/methods/mmddrift.html) on graph embeddings to flag **drift in the input data**.\n",
    "\n",
    "- [Learned Kernel detector](https://docs.seldon.io/projects/alibi-detect/en/stable/cd/methods/learnedkerneldrift.html) which flags **drift in the input data** using a (deep) learned kernel. The method trains a (deep) kernel on part of the data to maximise an estimate of the test power. Once the kernel is learned a permutation test is performed in the usual way on the value of the Maximum Mean Discrepancy (MMD) on the held out test set.\n",
    "\n",
    "- [Kolmogorov-Smirnov detector](https://docs.seldon.io/projects/alibi-detect/en/stable/cd/methods/ksdrift.html) to see if **drift** occurred **on graph level statistics** such as the number of nodes, edges and the average clustering coefficient.\n",
    "\n",
    "\n",
    "## Dataset\n",
    "\n",
    "We will train a classification model and detect drift on the *ogbg-molhiv* dataset. The dataset contains molecular graphs with both atom features (*atomic number-1, chirality, node degree, formal charge, number of H bonds, number of radical electrons, hybridization, aromatic?, in a ring?*) and bond level properties (*bond type (e.g. single or double), bond stereo code, conjugated?*). The goal is to predict whether a molecule inhibits HIV virus replication or not, so the task is binary classification. \n",
    "\n",
    "The dataset is split using the *scaffold splitting* procedure. This means that the molecules are split based on their 2D structural framework. Structurally different molecules are grouped into different subsets (train, validation, test) which could mean that there is drift between the splits.\n",
    "\n",
    "The dataset is retrieved from the [Open Graph Benchmark](https://ogb.stanford.edu/docs/graphprop/#ogbg-mol) dataset collection.\n",
    "\n",
    "## Dependencies\n",
    "\n",
    "Besides `alibi-detect`, this example notebook also uses [PyTorch Geometric](https://pytorch-geometric.readthedocs.io/en/stable/notes/installation.html) and [OGB](https://ogb.stanford.edu/docs/home/), both of which can be installed via pip/conda."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "910f9936",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import os\n",
    "import torch\n",
    "\n",
    "def set_seed(seed: int) -> None:\n",
    "    torch.manual_seed(seed)\n",
    "    torch.cuda.manual_seed(seed)\n",
    "    np.random.seed(seed)\n",
    "\n",
    "set_seed(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd497b63",
   "metadata": {},
   "source": [
    "## Load and analyze data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "6de4bc4a",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ogb.graphproppred import PygGraphPropPredDataset\n",
    "from torch_geometric.data import DataLoader"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e2eeaaa1",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "dataset_name = 'ogbg-molhiv'\n",
    "batch_size = 32\n",
    "dataset = PygGraphPropPredDataset(name=dataset_name)\n",
    "split_idx = dataset.get_idx_split()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "799fc764",
   "metadata": {},
   "source": [
    "We set some samples apart to serve as the reference data for our drift detectors. Note that the allowed format of the reference data is very flexible and can be `np.ndarray` or `List[Any]`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "545def71",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of reference instances: 1000\n",
      "Number of H0 instances: 500\n"
     ]
    }
   ],
   "source": [
    "n_ref = 1000\n",
    "n_h0 = 500\n",
    "\n",
    "idx_tr = split_idx['train']\n",
    "idx_sample = np.random.choice(idx_tr.numpy(), size=n_ref + n_h0, replace=False)\n",
    "idx_ref, idx_h0 = idx_sample[:n_ref], idx_sample[n_ref:]\n",
    "x_ref = [dataset[i] for i in idx_ref]\n",
    "x_h0 = [dataset[i] for i in idx_h0]\n",
    "idx_tr = torch.from_numpy(np.setdiff1d(idx_tr, idx_sample))\n",
    "print(f'Number of reference instances: {len(x_ref)}')\n",
    "print(f'Number of H0 instances: {len(x_h0)}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "03bcb089",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of train, val and test batches: 982, 129 and 129\n"
     ]
    }
   ],
   "source": [
    "dl_tr = DataLoader(dataset[idx_tr], batch_size=batch_size, shuffle=True)\n",
    "dl_val = DataLoader(dataset[split_idx['valid']], batch_size=batch_size, shuffle=False)\n",
    "dl_te = DataLoader(dataset[split_idx['test']], batch_size=batch_size, shuffle=False)\n",
    "print(f'Number of train, val and test batches: {len(dl_tr)}, {len(dl_val)} and {len(dl_te)}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "51e64302",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Dataset: PygGraphPropPredDataset(41127):\n",
      "=============================================================\n",
      "Number of graphs: 41127\n",
      "Number of node features: 9\n",
      "Number of edge features: 3\n",
      "Number of classes: 2\n",
      "\n",
      "Example: Data(edge_attr=[40, 3], edge_index=[2, 40], x=[19, 9], y=[1, 1])\n",
      "=============================================================\n",
      "Number of nodes: 19\n",
      "Number of edges: 40\n",
      "Average node degree: 2.11\n",
      "Contains isolated nodes: False\n",
      "Contains self-loops: False\n",
      "Is undirected: True\n"
     ]
    }
   ],
   "source": [
    "ds = dataset\n",
    "print()\n",
    "print(f'Dataset: {ds}:')\n",
    "print('=============================================================')\n",
    "print(f'Number of graphs: {len(ds)}')\n",
    "print(f'Number of node features: {ds.num_node_features}')\n",
    "print(f'Number of edge features: {ds.num_edge_features}')\n",
    "print(f'Number of classes: {ds.num_classes}')\n",
    "\n",
    "i = 0\n",
    "d = ds[i]\n",
    "\n",
    "print(f'\\nExample: {d}')\n",
    "print('=============================================================')\n",
    "\n",
    "print(f'Number of nodes: {d.num_nodes}')\n",
    "print(f'Number of edges: {d.num_edges}')\n",
    "print(f'Average node degree: {d.num_edges / d.num_nodes:.2f}')\n",
    "print(f'Contains isolated nodes: {d.contains_isolated_nodes()}')\n",
    "print(f'Contains self-loops: {d.contains_self_loops()}')\n",
    "print(f'Is undirected: {d.is_undirected()}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61092348",
   "metadata": {},
   "source": [
    "Let's plot some graph summary statistics such as the distribution of the node degrees, number of nodes and edges as well as the clustering coefficients:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "f914ef12",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import networkx as nx\n",
    "from networkx.algorithms.cluster import clustering\n",
    "from torch_geometric.utils import degree, to_networkx\n",
    "from tqdm import tqdm\n",
    "from typing import Tuple\n",
    "\n",
    "\n",
    "def degrees_and_clustering(loader: DataLoader) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:\n",
    "    degrees, c_coeff, num_nodes, num_edges = [], [], [], []\n",
    "    for data in tqdm(loader):\n",
    "        row, col = data.edge_index\n",
    "        deg = degree(row, data.x.size(0), dtype=data.x.dtype)\n",
    "        degrees.append(deg.numpy())\n",
    "        g = to_networkx(data, node_attrs=['x'], edge_attrs=['edge_attr'], to_undirected=True)\n",
    "        c = list(clustering(g).values())\n",
    "        c_coeff.append(c)\n",
    "        num_nodes += [d.num_nodes for d in data.to_data_list()]\n",
    "        num_edges += [d.num_edges for d in data.to_data_list()]\n",
    "    degrees = np.concatenate(degrees, axis=0)\n",
    "    c_coeff = np.concatenate(c_coeff, axis=0)\n",
    "    return degrees, c_coeff, np.array(num_nodes), np.array(num_edges)\n",
    "\n",
    "\n",
    "# x: nodes, edges, degree, cluster\n",
    "def plot_histogram(x: str, bins: int = None, log: bool = True) -> None:\n",
    "    if x == 'nodes':\n",
    "        vals = [num_nodes_tr, num_nodes_val, num_nodes_te]\n",
    "    elif x == 'edges':\n",
    "        vals = [num_edges_tr, num_edges_val, num_edges_te]\n",
    "    elif x == 'degree':\n",
    "        vals = [degree_tr, degree_val, degree_te]\n",
    "    elif x == 'cluster':\n",
    "        vals = [cluster_tr, cluster_val, cluster_te]\n",
    "    labels = ['train', 'val', 'test']\n",
    "    for v, l in zip(vals, labels):\n",
    "        plt.hist(v, density=True, log=log, label=l, bins=bins)\n",
    "    plt.title(f'{x} distribution')\n",
    "    plt.legend()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "2e5ff45c",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 982/982 [00:19<00:00, 50.56it/s]\n",
      "100%|██████████| 129/129 [00:03<00:00, 42.44it/s]\n",
      "100%|██████████| 129/129 [00:02<00:00, 48.92it/s]\n"
     ]
    }
   ],
   "source": [
    "degree_tr, cluster_tr, num_nodes_tr, num_edges_tr = degrees_and_clustering(dl_tr)\n",
    "degree_val, cluster_val, num_nodes_val, num_edges_val = degrees_and_clustering(dl_val)\n",
    "degree_te, cluster_te, num_nodes_te, num_edges_te = degrees_and_clustering(dl_te)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "d1d3b2ae",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average number and stdev of nodes, edges, degree and clustering coefficients:\n",
      "\n",
      "Train...\n",
      "Nodes: 25.3 +- 12.0\n",
      "Edges: 54.1 +- 26.1\n",
      "Degree: 2.1 +- 0.8\n",
      "Clustering: 0.001 +- 0.028\n",
      "\n",
      "Validation...\n",
      "Nodes: 27.8 +- 12.8\n",
      "Edges: 61.1 +- 27.8\n",
      "Degree: 2.2 +- 0.8\n",
      "Clustering: 0.003 +- 0.042\n",
      "\n",
      "Test...\n",
      "Nodes: 25.3 +- 12.5\n",
      "Edges: 55.6 +- 27.1\n",
      "Degree: 2.2 +- 0.7\n",
      "Clustering: 0.004 +- 0.058\n"
     ]
    }
   ],
   "source": [
    "print('Average number and stdev of nodes, edges, degree and clustering coefficients:')\n",
    "print('\\nTrain...')\n",
    "print(f'Nodes: {num_nodes_tr.mean():.1f} +- {num_nodes_tr.std():.1f}')\n",
    "print(f'Edges: {num_edges_tr.mean():.1f} +- {num_edges_tr.std():.1f}')\n",
    "print(f'Degree: {degree_tr.mean():.1f} +- {degree_tr.std():.1f}')\n",
    "print(f'Clustering: {cluster_tr.mean():.3f} +- {cluster_tr.std():.3f}')\n",
    "\n",
    "print('\\nValidation...')\n",
    "print(f'Nodes: {num_nodes_val.mean():.1f} +- {num_nodes_val.std():.1f}')\n",
    "print(f'Edges: {num_edges_val.mean():.1f} +- {num_edges_val.std():.1f}')\n",
    "print(f'Degree: {degree_val.mean():.1f} +- {degree_val.std():.1f}')\n",
    "print(f'Clustering: {cluster_val.mean():.3f} +- {cluster_val.std():.3f}')\n",
    "\n",
    "print('\\nTest...')\n",
    "print(f'Nodes: {num_nodes_te.mean():.1f} +- {num_nodes_te.std():.1f}')\n",
    "print(f'Edges: {num_edges_te.mean():.1f} +- {num_edges_te.std():.1f}')\n",
    "print(f'Degree: {degree_te.mean():.1f} +- {degree_te.std():.1f}')\n",
    "print(f'Clustering: {cluster_te.mean():.3f} +- {cluster_te.std():.3f}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "93b7ddf9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEICAYAAABcVE8dAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAZ7klEQVR4nO3dfZRU9Z3n8fdH7NCiiAg+8ZTuCa6isKK2iquTxc04gooYcxSfzpiNY8dRJw+7cdNOYgbPZFf27BnjeHxa3HDUieCwmqwSyBGdQXEmJCoeJsCAARVD4wMPCgEjRuN3/6jbWP1QdHX37arbdT+vc/pYde+te791LT7969/91e8qIjAzs9p3QLULMDOzynDgm5nlhAPfzCwnHPhmZjnhwDczywkHvplZTjjwbUCR9KykP+/nY8yW9KPk8ThJeyQNSmnf90u6NXk8VVJrGvtN9vfHkl5Ja39Wew6sdgFmWRYRvwEO6W47SV8G/jwizu5mf9enVBqSAjg2IjYm+34eOC6t/VvtcQvfrELS+ivBrLcc+NbvJG2S9C1Jv5K0S9I/SKovWn+dpI2S3pX0pKRRRevOlbQ+ed3dgDrs+yuS1kl6T9JTkj6bLJekH0jaKum3klZLmliivkZJz0naLelpYGTRugZJIenA5PmXJb2WbPu6pKskTQDuB85Mun92Jts+KOk+SUskvQ+ckyz7fofj/5Wk7cl5uqpoebvuq+TY/5w8Xp4s/tfkmLM6dhFJmpDsY6ektZIuKlr3oKR7JC1O3ssvJX2uu/+XNrA58K1SLgOmAY3Avwe+DCDpPwG3J+uPAd4AHk3WjQR+DHyXQgi/CpzVtkNJM4G/Ai4BjgCeBxYkq/8U+Dzw74Bhyf53lKhtPrAyOcbfANd0tZGkg4G7gOkRMRT4D8CqiFgHXA+siIhDIuKwopddCfx3YCjwz13s9ujkuKOT486V1G23TER8Pnl4UnLMf+hQax2wCFgKHAn8JfBIh31fDtwGDAc2JnVaDXPgW6XcFRFvRsS7FIJocrL8KmBeRLwcER8Ct1BoKTcA5wNrI+KxiPgIuBN4u2if1wO3R8S6iPgY+B/A5KSV/xGFkD0eULLNWx2LkjQOOA24NSI+jIjlSX2lfAJMlHRQRLwVEWu7ed9PRMS/RMQnEbG3xDZtx34OWEzhl1NfTaFw7WFORPw+Iv4J+ClwRdE2P4mIF5Jz9wif/j+xGuXAt0opDurf8emF0FEUWvUARMQeCi3x0cm6zUXrovg58Fng75Iui53AuxS6fEYnAXc3cA+wVdJcSYd2Udco4L2IeL9o2RtdbEeyzSwKv2jeSrpDju/mfW/uZn1Xxx5VauMeGAVsjohPOux7dNHzUv9PrEY58K3a3qQQ3MC+bpMRwBbgLWBs0ToVP6cQpl+NiMOKfg6KiJ8DRMRdEXEqcAKFrp2buzj+W8Dw5LhtxpUqNiKeiohzKXQ/rQceaFtV6iWl9pXo6thvJo/fB4YUrTu6m30VexMYK6n43/g4CufVcsqBb9W2APjPkiZLGkyhW+aXEbGJQvfGiZIuSS6afo32oXc/cIukEwEkDZN0afL4NElnJH3Z7wN7KXTHtBMRbwAvAbdJ+oyks4EZXRUq6ShJM5OA/hDYU7TPd4Axkj7Ti3PQduw/Bi4E/m+yfBVwiaQhksYD13Z43TvAH5XY5y8ptNr/m6Q6SVOT9/VoL+qzGuHAt6qKiGeAW4HHKbS2P0fhYiIRsR24FJhDoZvnWOBfil77E+B/Ao9K+i2wBpierD6UQuv7PQpdGTuA/1WijCuBMyh0Cf018HCJ7Q4A/guF1vO7wH8E/iJZ90/AWuBtSdvLff8UulXeS/b5CHB9RKxP1v0A+D2FYH8oWV9sNvBQ0qXVrt8/In5PIeCnA9uBe4E/K9q35ZB8AxQzs3xwC9/MLCcc+GZmOeHANzPLCQe+mVlOZHq2zJEjR0ZDQ0O1yzAzG1BWrly5PSKO6Lg804Hf0NDASy+9VO0yzMwGFEldfls8k106kmZImrtr165ql2JmVjMyGfgRsSgimocNG1btUszMakYmA9/MzNKXyT58STOAGePHj692KWY2wHz00Ue0trayd2+p2ahrR319PWPGjKGurq6s7TM9tUJTU1P4oq2Z9cTrr7/O0KFDGTFiBIUJVmtTRLBjxw52795NY2Nju3WSVkZEU8fXuEvHzGrK3r17az7sASQxYsSIHv0l48A3s5pT62HfpqfvM5OB72GZZmbpy+RF24hYBCxqamq6rtq1tJn00KROy1Zfs7oKlZhZTzS0LE51f5vmXNDtNjt37mT+/PnccMMNPdr3+eefz/z58znssMN6W95+ZbKFb2Y2kO3cuZN777230/KPP/54v69bsmRJv4U9ZLSFP6DM7vDlsNnuhjLLu5aWFl599VUmT55MXV0d9fX1DB8+nPXr1/PrX/+aiy++mM2bN7N3716+/vWv09zcDHw6ncyePXuYPn06Z599Nj//+c8ZPXo0TzzxBAcddFCf6spk4A+UcfgNLYvZVN9hoX8BmOXenDlzWLNmDatWreLZZ5/lggsuYM2aNfuGT86bN4/DDz+cDz74gNNOO40vfelLjBgxot0+NmzYwIIFC3jggQe47LLLePzxx7n66qv7VFcmu3QGytQKm+qvrHYJZjYAnH766e3Gyt91112cdNJJTJkyhc2bN7Nhw4ZOr2lsbGTy5MkAnHrqqWzatKnPdWSyhT9QTGoc12nZ6td/0+558QWjci72mFntOfjgg/c9fvbZZ3nmmWdYsWIFQ4YMYerUqV2OpR88ePC+x4MGDeKDDz7ocx0O/JR1/CUwlBZ2r5tTpWrMrBqGDh3K7t27u1y3a9cuhg8fzpAhQ1i/fj2/+MUvKlaXA78bbS30oROqXIiZ9Uo1/rIeMWIEZ511FhMnTuSggw7iqKOO2rdu2rRp3H///UyYMIHjjjuOKVOmVKyuTAb+QLloa2ZWyvz587tcPnjwYH72s591ua6tn37kyJGsWbNm3/JvfetbqdTki7ZmZjmRyRZ+rdk3mmc2HqZpZlWTyRa+mZmlz4FvZpYTDnwzs5zIZB9+lkbptPW/T6Lzl6zMzAaSTAZ+FqdHNrMBquP8Vn3eX/oDLw455BD27NmT+n47ymTgV1u7ue+7mD4hVZ5szcwqxIFvZpaylpYWxo4dy4033gjA7NmzOfDAA1m2bBnvvfceH330Ed///veZOXNmRety4CeKJznzNApm1hezZs3iG9/4xr7AX7hwIU899RRf+9rXOPTQQ9m+fTtTpkzhoosuquj9dx34ZmYpO/nkk9m6dStvvvkm27ZtY/jw4Rx99NF885vfZPny5RxwwAFs2bKFd955h6OPPrpidTnwK6DdDJrJ9QHfD9estl166aU89thjvP3228yaNYtHHnmEbdu2sXLlSurq6mhoaOhyWuT+5MCncJG20t04bV1Ine6YZWY1YdasWVx33XVs376d5557joULF3LkkUdSV1fHsmXLeOONNypekwPfzGpblUa+nXjiiezevZvRo0dzzDHHcNVVVzFjxgwmTZpEU1MTxx9/fMVrymTgZ+mLV/3Ft0c0q32rV3/adTty5EhWrFjR5XaVGIMPGQ38PHzxqpzbI5qZpSmTgZ9rxV/E8pewzCxFnjzNzCwnHPhmZjnhwDczywkHvplZTviirZnVtHaz36agnG/J79y5k/nz53PDDTf0eP933nknzc3NDBkypDfl7Zdb+GZmKdu5cyf33ntvr15755138rvf/S7ligrcwjczS1lLSwuvvvoqkydP5txzz+XII49k4cKFfPjhh3zxi1/ktttu4/333+eyyy6jtbWVP/zhD9x666288847vPnmm5xzzjmMHDmSZcuWpVqXA9/MLGVz5sxhzZo1rFq1iqVLl/LYY4/xwgsvEBFcdNFFLF++nG3btjFq1CgWLy7Mq7Vr1y6GDRvGHXfcwbJlyxg5cmTqdblLx8ysHy1dupSlS5dy8sknc8opp7B+/Xo2bNjApEmTePrpp/n2t7/N888/z7BhKd+KsQsVbeFLuhi4ADgU+GFELK3k8c3MKi0iuOWWW/jqV7/aad3LL7/MkiVL+O53v8sXvvAFvve97/VrLWW38CXNk7RV0poOy6dJekXSRkkt+9tHRPy/iLgOuB6Y1buSc2T2sPY/ZjYgDB06lN27dwNw3nnnMW/evH0TpG3ZsmXfzVGGDBnC1Vdfzc0338zLL7/c6bVp60kL/0HgbuDhtgWSBgH3AOcCrcCLkp4EBgG3d3j9VyJia/L4u8nrzMz6VTVuNjRixAjOOussJk6cyPTp07nyyis588wzATjkkEP40Y9+xMaNG7n55ps54IADqKur47777gOgubmZadOmMWrUqNQv2ioiyt9YagB+GhETk+dnArMj4rzk+S0AEdEx7NteL2AO8HREPFNim2agGWDcuHGnVuImAWmP0+2tbmfL9GRqZt1at24dEybk58bUXb1fSSsjoqnjtn3twx8NbC563gqcsZ/t/xL4E2CYpPERcX/HDSJiLjAXoKmpqfzfRnnQsVvHvwDMrAcqetE2Iu4C7upuuzzcAMXMrNL6GvhbgLFFz8cky/okDzdAKZdvlGLWcxFBoQe5tvWkSx76HvgvAsdKaqQQ9JcDvndfL3UV7mbWM/X19ezYsYMRI0bUdOhHBDt27KC+vr7s15Qd+JIWAFOBkZJagb+OiB9Kugl4isLInHkRsbZnZXd5LHfpmFmvjBkzhtbWVrZt21btUvpdfX09Y8aMKXv7sgM/Iq4osXwJsKTsI5Z3LHfpmFmv1NXV0djYWO0yMslTK5iZ5UQmA1/SDElzd+3ysEMzs7RkcrZMd+mUyePyzawHMtnCNzOz9GUy8N2lY2aWvkwGfkQsiojmSswPbWaWF5nsw7f987dvzaw3MtnCNzOz9GUy8N2Hb2aWvkwGvvvwzczS5z78GjGpcRwU3cilGnf5MbNsy2QL38zM0ufANzPLiUx26fTn9MhZuX9tv5s9zFMtmFk7mQx8z6WTEs+1Y2ZF3KVjZpYTDnwzs5xw4JuZ5YQD38wsJzJ50db6rssJ1qpQh5llRyZb+J5Lx8wsfZls4XtYZuV09b0ET8tgVpsy2cI3M7P0OfDNzHLCgW9mlhMOfDOznHDg51nHuXbMrKY58M3McsKBb2aWE5kMfH/xyswsfZkMfN/EvMpmD3P/vlkNyuQ3ba269s3D45uim9WUTLbwzcwsfW7h50hDy+J2zzfVV6kQM6sKt/DNzHLCLfwcGTqhpd3zSXSeM7+k2cN8E3SzAc4tfDOznHDgm5nlhLt0rCyTGse1G6YJHqppNtC4hW9mlhNu4Vvvdfw2ri/qmmVaxVr4kiZIul/SY5L+olLHNTOzgrICX9I8SVslremwfJqkVyRtlNRS6vUAEbEuIq4HLgPO6n3JZmbWG+W28B8EphUvkDQIuAeYDpwAXCHpBEmTJP20w8+RyWsuAhYDS1J7B2ZmVpay+vAjYrmkhg6LTwc2RsRrAJIeBWZGxO3AhSX28yTwpKTFwPyutpHUDDQDjBvXgy8GmZnZfvXlou1oYHPR81bgjFIbS5oKXAIMZj8t/IiYC8wFaGpqij7UZ2ZmRSo2SicingWeLWdbSTOAGePHj+/PkszMcqUvo3S2AGOLno9JlvWZb4BiZpa+vrTwXwSOldRIIegvB65MpSobmDwu3yzTyh2WuQBYARwnqVXStRHxMXAT8BSwDlgYEWvTKMr3tDUzS1+5o3SuKLF8Cf0wxDIiFgGLmpqarkt732ZmeeW5dMzMciKTge8uHTOz9GVy8jR36eSEL/KaVVQmW/hmZpa+TAa+u3TMzNLnLh3rP+6yMcuUTAa+DQyTGsub3G7167/p50rMrByZ7NIxM7P0ZTLw3YdvZpa+TAa+J08zM0tfJgPfzMzS54u2VjENLYv3Pd4054IqVmKWT27hm5nlRCYD3xdtzczSl8nA90VbM7P0ZTLwrfYV9+fvb5mZpceBb2aWEx6lYxWzqb6PtzzuODdPp/W+5mO2P27hm5nlhAPfzCwnMhn4HpZpZpa+TAa+h2WamaXPF22t35U7bz7rOrzuoUntF5TYj+fbNytPzQZ+qTHdQydUuBAr26b6K2F20YJyf1GU4lE9Zu1kskvHzMzS58A3M8sJB76ZWU448M3McsKBb2aWE5kMfH/xyswsfZkMfH/xyswsfZkMfDMzS1/NfvHKrMc6flGruy9m+YtdNsC4hW9mlhMOfDOznHDgm5nlhAPfzCwnfNHWatKkhyZ1P9tmd9MvF61ffc3q7o/ZxevLeZ1ZpbiFb2aWEw58M7OccOCbmeVERQNf0sGSXpJ0YSWPa2ZmZV60lTQPuBDYGhETi5ZPA/4OGAT8n4iY082uvg0s7GWtZvtVfFtL38rSrLNyR+k8CNwNPNy2QNIg4B7gXKAVeFHSkxTC//YOr/8KcBLwb0B930o2M7PeKCvwI2K5pIYOi08HNkbEawCSHgVmRsTtFP4aaEfSVOBg4ATgA0lLIuKTLrZrBpoBxo3r402szcxsn76Mwx8NbC563gqcUWrjiPgOgKQvA9u7Cvtku7nAXICmpqboQ31mZlak4l+8iogHu9tG0gxgxvjx4/u/IDOznOjLKJ0twNii52OSZX3mG6CYmaWvLy38F4FjJTVSCPrLgStTqcqsFzbVf/rxm0RGrv8Uz5nf0/n1+zqffnfz9Xfa3vP317qyWviSFgArgOMktUq6NiI+Bm4CngLWAQsjYm0aRfmetmZm6St3lM4VJZYvAZakWlFhv4uARU1NTdelvW8zs7zy1ApmZjmRyemRPUonnzpNL9zPr+vRMTpOpQyUM/Fxu9qSfZQ7ZXKXx8z5dMs+J32TyRa+R+mYmaUvk4FvZmbpy2Tge5SOmVn6Mhn47tIxM0tfJgPfzMzS58A3M8uJTAa++/DNzNKXycB3H76ZWfoyGfhmZpY+B76ZWU448M3MciKTge+LtmZm6ctk4PuirZlZ+jI5W6ZZLWtoWdzu+aY5F/TPfutT2a3VkEy28M3MLH0OfDOznHDgm5nlhAPfzCwnMnnR1rc4NCut48XZaitVT1oXoy09mWzhe1immVn6Mhn4ZmaWPge+mVlOOPDNzHLCgW9mlhMOfDOznMjksEyzPGkb1ui5b6y/ZbKF7+mRzczSl8nA9zh8M7P0ZTLwzcwsfQ58M7OccOCbmeWEA9/MLCcc+GZmOeHANzPLCX/xyszKkrV5+LMo6/cGcAvfzCwnHPhmZjnhwDczywkHvplZTlQs8CVNlfS8pPslTa3Ucc3MrKCswJc0T9JWSWs6LJ8m6RVJGyW1dLObAPYA9UBr78o1M7PeKndY5oPA3cDDbQskDQLuAc6lEOAvSnoSGATc3uH1XwGej4jnJB0F3AFc1bfSzcysJ8oK/IhYLqmhw+LTgY0R8RqApEeBmRFxO3Dhfnb3HjC41EpJzUAzwLhx48opz8zMytCXPvzRwOai563Jsi5JukTS/wb+nsJfC12KiLkR0RQRTUcccUQfyjMzs2IV+6ZtRPwY+HE520qaAcwYP358/xZlZpYjfWnhbwHGFj0fkyzrM9/xyswsfX0J/BeBYyU1SvoMcDnwZDplmZlZ2sodlrkAWAEcJ6lV0rUR8TFwE/AUsA5YGBFr0yjKNzE3M0tfuaN0riixfAmwJNWKCvtdBCxqamq6Lu19m5nlladHNhtA+nOK4rZ99/dUvvt7D6WO3faaoRP6dozu9l+unp6j/t5/uTI5l467dMzM0pfJwPcoHTOz9GUy8M3MLH2ZDHx36ZiZpS+Tge8uHTOz9GUy8M3MLH0OfDOznMhk4LsP38wsfZkMfPfhm5mlTxFR7RpKkrQNeKMHLxkJbO+ncgYyn5fOfE468znpbKCek89GRKcbimQ68HtK0ksR0VTtOrLG56Uzn5POfE46q7VzkskuHTMzS58D38wsJ2ot8OdWu4CM8nnpzOekM5+TzmrqnNRUH76ZmZVWay18MzMrwYFvZpYTNRP4kqZJekXSRkkt1a6nWiRtkrRa0ipJLyXLDpf0tKQNyX+HV7vO/iRpnqStktYULevyHKjgruRz8ytJp1Sv8v5V4rzMlrQl+bysknR+0bpbkvPyiqTzqlN1/5E0VtIySf8maa2kryfLa/azUhOBL2kQcA8wHTgBuELSCdWtqqrOiYjJReOHW4B/jIhjgX9MnteyB4FpHZaVOgfTgWOTn2bgvgrVWA0P0vm8APwg+bxMTu5TTfLv53LgxOQ19yb/zmrJx8B/jYgTgCnAjcn7rtnPSk0EPnA6sDEiXouI3wOPAjOrXFOWzAQeSh4/BFxcxVr6XUQsB97tsLjUOZgJPBwFvwAOk3RMZSqtrBLnpZSZwKMR8WFEvA5spPDvrGZExFsR8XLyeDewDhhNDX9WaiXwRwObi563JsvyKIClklZKak6WHRURbyWP3waOqk5pVVXqHPizAzclXRTzirr7cnVeJDUAJwO/pIY/K7US+PapsyPiFAp/ft4o6fPFK6MwDjfXY3F9Dtq5D/gcMBl4C/jb6pZTeZIOAR4HvhERvy1eV2uflVoJ/C3A2KLnY5JluRMRW5L/bgV+QuHP8Hfa/vRM/ru1ehVWTalzkOvPTkS8ExF/iIhPgAf4tNsmF+dFUh2FsH8kIn6cLK7Zz0qtBP6LwLGSGiV9hsLFpierXFPFSTpY0tC2x8CfAmsonItrks2uAZ6oToVVVeocPAn8WTICYwqwq+jP+ZrXoQ/6ixQ+L1A4L5dLGiypkcKFyhcqXV9/kiTgh8C6iLijaFXtflYioiZ+gPOBXwOvAt+pdj1VOgd/BPxr8rO27TwAIyiMNtgAPAMcXu1a+/k8LKDQPfERhX7Wa0udA0AURni9CqwGmqpdf4XPy98n7/tXFALtmKLtv5Ocl1eA6dWuvx/Ox9kUumt+BaxKfs6v5c+Kp1YwM8uJWunSMTOzbjjwzcxywoFvZpYTDnwzs5xw4JuZ5YQD38wsJxz4ZmY58f8Bgqv7ZGan+2EAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEICAYAAABcVE8dAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAalUlEQVR4nO3dfZQV1Znv8e9PbGmRvoDgG29pIi6DsRM0HV+WrnVJHCMakSST+EZmyOi1zVXnmizN2M6owbsyI3Mzy1FWNIYYLiaKyviSgOKIJijOkkSBkAECCmorDSpIhKARB8xz/zjV3MOhG5ru6nOqT/0+a/Xi1K5dVc9umofdu3btUkRgZmbV74BKB2BmZuXhhG9mlhNO+GZmOeGEb2aWE074ZmY54YRvZpYTTviWKZLGSWqtcAzfkPQfRdvvSfp4Suf+e0l3J5/rJYWkA1M698gk1j5pnM+qjxO+2T5ERP+IeHVvdTr7H1VE/FNE/I804pLUIukvis79RhLrR2mc36qPE75ZmaTVkzfrKid863GShkp6WNImSa9J+l9F+w6WNFPSu5J+D3y25NgTJf1W0jZJ/ybpQUnfK9p/rqRlkrZIel7Sp4r2XSdpfXLsS5LO6CC+wZLmSPqjpBeAo0v2h6TRyedzJP0+Oed6SddKOgR4AhiaDKm8l7R5iqSHJN0r6Y/AN5Kye0tCuETSBklvSrq26LozS9q667cIST8DRgJzk+v9XekQURLDHEl/kLRW0mVF55oiabaknyZtWSmpce9/k9bbOeFbj5J0ADAX+B0wDDgD+Jaks5Iq36WQYI8GzgImFx17EPAoMBM4FLgf+HLR/hOAGcDlwGDgR8AcSX0lHQtcBXw2IuqSc7d0EOYdwHbgKOCS5KsjPwEuT855PPCriHgfOBvYkAyp9I+IDUn9icBDwEDgvg7O+TngGOALwHXFwzQdiYi/At4AJiTX+z/tVHsAaAWGAl8F/knS54v2n5fUGQjMAX6wr+ta7+aEbz3ts8BhEfG/I+K/krHwHwMXJvvPB/4xIv4QEeuAaUXHngIcCEyLiB0R8QjwQtH+JuBHEfGbiPgoIu4BPkyO+wjoCxwnqSYiWiLildLgkhucfwncFBHvR8QK4J69tGdHcs7/FhHvRsTSfbR/UUT8PCL+HBEfdFDn5uTay4H/C1y0j3Puk6QRwGnAdRGxPSKWAXcDf11U7T8iYl4y5v8z4NPdva5lmxO+9bSPURjq2NL2Bfw9cESyfyiwrqj+60WfhwLrY/cV/orrfgy4puTcI4ChEbEW+BYwBdgo6QFJQ9uJ7zAK/6l0FEOpvwTOAV6X9KykU/dStzTeztR5nUK7u2so8IeI2FZy7mFF228Vff4TUOv7DNXNCd962jrgtYgYWPRVFxHnJPvfpJCk24ws+vwmMEySisqK666j8NtB8bn7RcT9ABExKyJOp/AfQwD/3E58m4Cde4lhNxHxYkRMBA4Hfg7MbtvV0SEdnauDNo0E2oaD3gf6Fe07cj/OvQE4VFJdybnXdyIeq1JO+NbTXgC2JTdQD5bUR9Lxktpuzs4Grpc0SNJw4G+Ljl1EYWjmKkkHSpoInFS0/8fANyWdrIJDJH1RUp2kYyV9XlJfCuPzHwB/Lg0uGc54BJgiqZ+k4yi6j1BM0kGSJkkaEBE7gD8WnfNtYLCkAV34Ht2YXPuTwN8ADybly4BzJB0q6UgKv7EUexto9/mAZHjseeAWSbXJzexLgdIbxpYjTvjWo5KEei4wFngNeIfCWHJbYryZwlDDa8B8CmPJbcf+F/AVColqC/B14DEK4/RExGLgMgo3G98F1gLfSA7vC0xNrvcWhR759R2EeRXQP6k3k8I4ekf+CmhJZt18E5iUxLKawk3lV5Phpf0Zlnk2if2XwL9ExPyk/GcUbna3UPjePFhy3C3ADcn1rmVPFwH1FHr7jwLfjYin9yMuqzLyC1CsN5H0G+CuiNhbUjazdriHb5km6b9LOjIZ0pkMfAr490rHZdYb+Y68Zd2xFMb5DwFeBb4aEW9WNiSz3slDOmZmOeEhHTOznMj0kM6QIUOivr6+0mGYmfUqS5YseSciDistz2TClzQBmDB69GgWL15c6XDMzHoVSe0+LZ7JIZ2ImBsRTQMGdOUZFjMza08mE76ZmaXPCd/MLCcyOYZvZtZVO3bsoLW1le3bt1c6lB5XW1vL8OHDqamp6VR9J3wzqyqtra3U1dVRX1/P7gutVpeIYPPmzbS2tjJq1KhOHZPJIR1JEyRN37p1a6VDMbNeZvv27QwePLiqkz2AJAYPHrxfv8lkMuF7lo6ZdUe1J/s2+9vOTCZ8MzNLn8fwu6jhnoY9ypZPXl6BSMxsb+qbH0/1fC1Tv7jPOlu2bGHWrFlcccUV+3Xuc845h1mzZjFw4MCuhrdX7uGbmaVsy5Yt3HnnnXuU79y5c6/HzZs3r8eSPbiHn6q2nkRnegBmVr2am5t55ZVXGDt2LDU1NdTW1jJo0CBWr17Nyy+/zJe+9CXWrVvH9u3bufrqq2lqagKgvr6exYsX895773H22Wdz+umn8/zzzzNs2DB+8YtfcPDBB3crLvfwU9RSezEttRfDFN9sNsuzqVOncvTRR7Ns2TK+//3vs3TpUm6//XZefvllAGbMmMGSJUtYvHgx06ZNY/PmzXucY82aNVx55ZWsXLmSgQMH8vDDD3c7Lvfwzcx62EknnbTbXPlp06bx6KOPArBu3TrWrFnD4MGDdztm1KhRjB07FoDPfOYztLS0dDsOJ/xOau8mrZlZZxxyyCG7Pj/zzDM8/fTTLFq0iH79+jFu3Lh259L37dt31+c+ffrwwQcfdDuOTA7p+MErM+vN6urq2LZtW7v7tm7dyqBBg+jXrx+rV6/m17/+ddniymQPPyLmAnMbGxsvq3QsbTdi68ZUOBAz65JKTKIYPHgwp512GscffzwHH3wwRxxxxK5948eP56677mLMmDEce+yxnHLKKWWLK5MJvxq0N/fXs3fM8mPWrFntlvft25cnnnii3X1t4/RDhgxhxYoVu8qvvfbaVGJywk9Rw6iRuz7X0cy2VVMrGI2Z2e6c8HtQ3Zjm3bYb7mn207hmVjGZvGlrZmbpc8IvtykD/GCWmVWEE76ZWU444ZuZ5YRv2lZK6bDOFD9kZtYj0h5C7YF/q/379+e9995L/bylnPD3oaX2YgAaGLmPmmZm2eaEb2aWsubmZkaMGMGVV14JwJQpUzjwwANZsGAB7777Ljt27OB73/seEydOLGtcZR3Dl/QlST+W9KCkL5Tz2mZm5XLBBRcwe/bsXduzZ89m8uTJPProoyxdupQFCxZwzTXXEBFljavTPXxJM4BzgY0RcXxR+XjgdqAPcHdEdPh4aUT8HPi5pEHAvwDzuxp4b1X8NG6b5a+9UYFIzKynnHDCCWzcuJENGzawadMmBg0axJFHHsm3v/1tFi5cyAEHHMD69et5++23OfLII8sW1/4M6cwEfgD8tK1AUh/gDuBMoBV4UdIcCsn/lpLjL4mIjcnnG5LjzMyq0te+9jUeeugh3nrrLS644ALuu+8+Nm3axJIlS6ipqaG+vr7dZZF7UqcTfkQslFRfUnwSsDYiXgWQ9AAwMSJuofDbwG4kCZgKPBERS9u7jqQmoAlg5MjK3Cjdbe37dnrkZmb7csEFF3DZZZfxzjvv8OyzzzJ79mwOP/xwampqWLBgAa+//nrZY+ruTdthwLqi7Vbg5L3U/1vgL4ABkkZHxF2lFSJiOjAdoLGxsbwDXGZWfSo05fmTn/wk27ZtY9iwYRx11FFMmjSJCRMm0NDQQGNjI5/4xCfKHlNZZ+lExDRgWjmv2VnFyxl77XszS8Py5f9/scQhQ4awaNGiduuVYw4+dH+WznpgRNH28KSsW/zGKzOz9HU34b8IHCNplKSDgAuBOd0NKiLmRkTTgAE5WmSsbVE1L6xmZj2k0wlf0v3AIuBYSa2SLo2IncBVwJPAKmB2RKzsblDu4ZuZpW9/Zulc1EH5PGBeahFR/nfaNtzT4HF7M6t6Xi3TzCwnMrmWjqQJwITRo0dXOpTK8EqaZtYDMpnwyz2kY2bVa7cHKVPQmfdSb9myhVmzZnHFFVfs9/lvu+02mpqa6NevX1fC2ysP6ZiZpWzLli3ceeedXTr2tttu409/+lPKERVksoef+yEdM+vVmpubeeWVVxg7dixnnnkmhx9+OLNnz+bDDz/ky1/+MjfffDPvv/8+559/Pq2trXz00UfceOONvP3222zYsIHPfe5zDBkyhAULFqQaVyYTft6GdEpX0PTqmWa929SpU1mxYgXLli1j/vz5PPTQQ7zwwgtEBOeddx4LFy5k06ZNDB06lMcfLzzlv3XrVgYMGMCtt97KggULGDJkSOpxeUjHzKwHzZ8/n/nz53PCCSdw4oknsnr1atasWUNDQwNPPfUU1113Hc899xzleNA0kz18K+FZO2a9VkRw/fXXc/nll++xb+nSpcybN48bbriBM844g5tuuqlHY8lkws/7GL5fkmLWu9XV1bFt2zYAzjrrLG688UYmTZpE//79Wb9+PTU1NezcuZNDDz2Ur3/96wwcOJC77757t2N7Ykgnkwk/b2P4ZtZzOjONMm2DBw/mtNNO4/jjj+fss8/m4osv5tRTTwWgf//+3Hvvvaxdu5bvfOc7HHDAAdTU1PDDH/4QgKamJsaPH8/QoUNTv2mrcr9TcX80NjbG4sWLe/w6ac/T7Qm79fA9pGPWoVWrVjFmTH7WSmmvvZKWRERjaV3ftDUzywknfDOznMjkGH7eb9ruk2ftmO1VRFB4hXZ1298h+Uz28HP5AhQzS0VtbS2bN2/e72TY20QEmzdvpra2ttPHZLKHb2bWVcOHD6e1tZVNmzZVOpQeV1tby/Dhwztd3wnfzKpKTU0No0aNqnQYmeSE30t4vR0z665MjuGbmVn6nPDNzHIikwlf0gRJ07du9XRDM7O0ZDLhe1qmmVn6MpnwzcwsfU74ZmY54YRvZpYTTvhmZjnhB696qd0exErW81/+2hteSM3MOuQevplZTjjhly41bGZWpTI5pNOT6+HXNz++23ZL51cWNTPr1TKZ8P0S865pGDVy13g+VOblzWaWXR7SMTPLCSd8M7OccMI3M8uJTI7hW0ram4HkefpmueUevplZTjjhm5nlhId0qljpe3ABlpcO83iIxyw3nPBzrqFo3n4bz983q065S/h1Y5p3225gz16wmVk18hi+mVlOlC3hSxoj6S5JD0n6n+W6rpmZFXQq4UuaIWmjpBUl5eMlvSRpraTmjo4HiIhVEfFN4HzgtK6HbGZmXdHZHv5MYHxxgaQ+wB3A2cBxwEWSjpPUIOmxkq/Dk2POAx4H5qXWAjMz65RO3bSNiIWS6kuKTwLWRsSrAJIeACZGxC3AuR2cZw4wR9LjwKz26khqApoARo70DdW0tTdVcw9tUzc9ZdOsqnRnls4wYF3RditwckeVJY0DvgL0ZS89/IiYDkwHaGxsjG7EZ2ZmRco2LTMingGe6UzdnnwBiu3brt8C2t6V63n5ZlWhO7N01gMjiraHJ2XdFhFzI6JpwAC/ftDMLC3dSfgvAsdIGiXpIOBCYE46YZmZWdo6Oy3zfmARcKykVkmXRsRO4CrgSWAVMDsiVqYRlKQJkqZv3eqbhmZmaensLJ2LOiifRw9MsfQ7bc3M0uelFczMciKTCd9DOmZm6ctkwvcsHTOz9OVueWRLUfHLVPxUrlnmZbKH7yEdM7P0ZbKH71k62eK3YplVh0z28M3MLH1O+GZmOZHJhO8xfDOz9GUy4XtapplZ+jKZ8M3MLH2ZnKVjVWBKyW9nnqdvVnHu4ZuZ5UQme/h+41X2NdzTAEXvx/WsfLPsy2QP3zdtzczSl8mEb2Zm6XPCNzPLCSd8M7OccMI3M8uJTCZ8L61gZpa+TCZ8z9IxM0tfJhO+mZmlzwnfzCwnnPDNzHIik0srWC9Uulhad+p6oTWzHuGEb2VR3/z4rs8ttRUMxCzHPKRjZpYTTvhmZjmRyYTvB6/MzNKXyYTvB6/MzNKXyYRvZmbpc8I3M8sJT8u0VDQUve6wPXU0s23V1DJFY2btccK37Cl9MGtvD2LtT900jzXrhTykY2aWE074ZmY54YRvZpYTHsO3TGn35u89DXsULZ+8vP3jiuqW1jHLO/fwzcxywj18K5uW2ovTO9n+LMecxjU8g8eqQFl7+JIOkbRY0rnlvK6ZmXUy4UuaIWmjpBUl5eMlvSRpraTmTpzqOmB2VwI1M7Pu6eyQzkzgB8BP2wok9QHuAM4EWoEXJc0B+gC3lBx/CfBp4PeAX39hZlYBnUr4EbFQUn1J8UnA2oh4FUDSA8DEiLgF2GPIRtI44BDgOOADSfMi4s/t1GsCmgBGjtz74/pmZtZ53blpOwxYV7TdCpzcUeWI+AcASd8A3mkv2Sf1pgPTARobG6OrwRW/Uq9Y3ZiuntHMrHcr+yydiJhZ7muamVn3ZumsB0YUbQ9PyrrNb7wyM0tfdxL+i8AxkkZJOgi4EJiTRlB+45WZWfo6NaQj6X5gHDBEUivw3Yj4iaSrgCcpzMyZEREr0whK0gRgwujRo9M4neVZOR7QKicv6Wzd0NlZOhd1UD4PmJdqRIXzzgXmNjY2Xpb2uc3M8spr6ZiZ5UQmE75v2pqZpS+Ti6d5SKc67eu9t6lf754G2Mc1l7/2RpmiMau8TPbwzcwsfU74ZmY5kcmE7zF8M7P0ZTLh+8ErM7P0ZTLhm5lZ+pzwzcxyIpMJ32P4Zmbpy2TC9xi+mVn6MpnwzcwsfU74ZmY5kcmlFcwyJ81liXtqiWMvnWz7kMkevm/ampmlL5MJ3zdtzczSl8mEb2Zm6XPCNzPLCd+0tV4prbX1dzvPPQ2770z2la6Z3zBq5J5190cH561mDe18v5ZPXl6BSPLNPXwzs5zIZML3LB0zs/RlMuF7lo6ZWfoymfDNzCx9TvhmZjnhhG9mlhNO+GZmOeGEb2aWE074ZmY54YRvZpYTmUz4fvDKzCx9mUz4fvDKzCx9mUz4ZmaWPid8M7OccMI3M8sJr4dvllH1zY/vUdZS27PnB2iZ+sX0LmKZ4h6+mVlOOOGbmeWEE76ZWU444ZuZ5YQTvplZTjjhm5nlRNkSvqRxkp6TdJekceW6rpmZFXQq4UuaIWmjpBUl5eMlvSRpraTmfZwmgPeAWqC1a+GamVlXdfbBq5nAD4CfthVI6gPcAZxJIYG/KGkO0Ae4peT4S4DnIuJZSUcAtwKTuhe6mZntj04l/IhYKKm+pPgkYG1EvAog6QFgYkTcApy7l9O9C/TtaKekJqAJYOTIkZ0Jz8zMOqE7SysMA9YVbbcCJ3dUWdJXgLOAgRR+W2hXREwHpgM0NjZGN+Iz6xU6WuIgrfN6qQRrU7a1dCLiEeCRztSVNAGYMHr06J4NyswsR7ozS2c9MKJoe3hS1m1+AYqZWfq6k/BfBI6RNErSQcCFwJx0wjIzs7R1dlrm/cAi4FhJrZIujYidwFXAk8AqYHZErEwjKL/T1swsfZ2dpXNRB+XzgHmpRlQ471xgbmNj42Vpn9vMLK+8tIKZWU5kMuF7SMfMLH2ZTPiepWNmlj6/09asyvXUg13d1e47e3P2kFi53yucyR6+h3TMzNKXyYTvIR0zs/RlMuGbmVn6nPDNzHIikwnfY/hmZunLZML3GL6ZWfoymfDNzCx9TvhmZjnhhG9mlhOZTPi+aWtmlr5MJnzftDUzS18mE76ZmaXPCd/MLCec8M3McsLLI5v1YsXL67bUlu9axdJcyjetpZw7iqmn21Du5Y73l3v4ZmY5kcmE72mZZmbpy2TC97RMM7P0ZTLhm5lZ+pzwzcxywgnfzCwnnPDNzHLCCd/MLCec8M3MckIRUekYOiRpE/B6Fw4dAryTcji9gdudH3lsM7jdnfWxiDistDDTCb+rJC2OiMZKx1Fubnd+5LHN4HZ39zwe0jEzywknfDOznKjWhD+90gFUiNudH3lsM7jd3VKVY/hmZranau3hm5lZCSd8M7OcqKqEL2m8pJckrZXUXOl40iRphqSNklYUlR0q6SlJa5I/ByXlkjQt+T78p6QTKxd590gaIWmBpN9LWinp6qS8qtsuqVbSC5J+l7T75qR8lKTfJO17UNJBSXnfZHttsr++kvF3h6Q+kn4r6bFku+rbDCCpRdJyScskLU7KUv05r5qEL6kPcAdwNnAccJGk4yobVapmAuNLypqBX0bEMcAvk20ofA+OSb6agB+WKcaesBO4JiKOA04Brkz+Xqu97R8Cn4+ITwNjgfGSTgH+GfjXiBgNvAtcmtS/FHg3Kf/XpF5vdTWwqmg7D21u87mIGFs05z7dn/OIqIov4FTgyaLt64HrKx1Xym2sB1YUbb8EHJV8Pgp4Kfn8I+Ci9ur19i/gF8CZeWo70A9YCpxM4WnLA5PyXT/zwJPAqcnnA5N6qnTsXWjr8CSxfR54DFC1t7mo7S3AkJKyVH/Oq6aHDwwD1hVttyZl1eyIiHgz+fwWcETyuSq/F8mv7CcAvyEHbU+GNpYBG4GngFeALRGxM6lS3LZd7U72bwUGlzfiVNwG/B3w52R7MNXf5jYBzJe0RFJTUpbqz/mBaUVqlRURIalq59hK6g88DHwrIv4oade+am17RHwEjJU0EHgU+ESFQ+pRks4FNkbEEknjKh1PBZweEeslHQ48JWl18c40fs6rqYe/HhhRtD08Katmb0s6CiD5c2NSXlXfC0k1FJL9fRHxSFKci7YDRMQWYAGF4YyBkto6asVt29XuZP8AYHOZQ+2u04DzJLUAD1AY1rmd6m7zLhGxPvlzI4X/4E8i5Z/zakr4LwLHJHf0DwIuBOZUOKaeNgeYnHyeTGF8u638r5M7+acAW4t+LexVVOjK/wRYFRG3Fu2q6rZLOizp2SPpYAr3LVZRSPxfTaqVtrvt+/FV4FeRDO72FhFxfUQMj4h6Cv9+fxURk6jiNreRdIikurbPwBeAFaT9c17pGxUp3/Q4B3iZwljnP1Q6npTbdj/wJrCDwnjdpRTGK38JrAGeBg5N6orCjKVXgOVAY6Xj70a7T6cwtvmfwLLk65xqbzvwKeC3SbtXADcl5R8HXgDWAv8G9E3Ka5Pttcn+j1e6Dd1s/zjgsby0OWnj75KvlW35K+2fcy+tYGaWE9U0pGNmZnvhhG9mlhNO+GZmOeGEb2aWE074ZmY54YRvZpYTTvhmZjnx/wBeq9sSlfS/kwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEICAYAAABcVE8dAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAZI0lEQVR4nO3df5BU5Z3v8fdHQEZ+rBhQkV8ZNlgGkBV0JHgp7zVrjGhAzC/xV9asLug1bjYpdR1zTUJu5S7WzZZlLDEuG1k0imYWNcE4rtxkUXZrSeRHSAFBAyKGAZWBBIJGjZjv/aMb0rYz0DN9uvtMn8+riqL76XOe821+fOaZ5zzztCICMzOrf0fVugAzM6sOB76ZWUY48M3MMsKBb2aWEQ58M7OMcOCbmWWEA9+qQtIiSd+sdR1dJenzkv6z4Pnrkv48ob6/Ium7+ceNkkJS74T6HpWvtVcS/Vl9cOCbdUFEDIiIrYc7RtI5ktpK6OsfIuJvkqhL0jZJHyvo+9f5Wt9Non+rDw5869GSGhFXW0+t23o2B75VhKRJktZK2i/p+0BD0evTJa2TtFfSf0n6i4LXTpf08/y5/yrp+wengw6OniXdIulV4F9K6G+YpEcltUt6SdIXD1P3YElLJf1O0nPAh4peD0lj8o8vlPTLfJ07JN0kqT/wFDAsP6Xyev76cyUtkfSgpN8Bn8+3PVhUwtWSdkp6RdJNBdd9z5RY4XcRkr4HjAKeyF/v74uniPI1LJX0G0lbJM0u6GuupBZJD+Tfy0ZJTYf567UeyoFviZN0NPAD4HvAB4B/BT5d8PokYCFwLTAY+CdgqaS++XMfBxblz30Y+GTRJYbmX/sgMOcI/R0FPAH8AhgOnAt8SdL5nZQ/H3gLOAm4Ov+rM/cB10bEQOBU4N8j4g3gAmBnfkplQETszB8/E1gCDAIe6qTPjwInAx8HbimcpulMRHwO+DUwI3+9/9vBYY8AbcAw4DPAP0j6y4LXL8ofMwhYCtx9pOtaz+PAt0qYAvQB7oyIdyJiCbCq4PU5wD9FxM8i4t2IuB94O3/eFKA3cFf+3MeA54r6/yPw9Yh4OyLePEJ/ZwLHR8T/jog/5Off/xm4tLjo/A3OTwNfi4g3ImIDcP9h3uc7wDhJfxYRv42ItUf4c1kZET+IiD/m6+7IN/LXXk/uu5fLjtDnEUkaCUwFbomItyJiHfBd4K8KDvvPiGjNz/l/Dzit3Ota+jjwrRKGATvivTvzvVzw+IPAjfnpl72S9gIj8+d1dO72ov7bI+KtEvv7ILnplcLXvgKc2EHdx5P7YlN4vZc7OO6gTwMXAi9LelbSWYc5tqP3caRjXib3Hso1DPhNROwv6nt4wfNXCx7/HmjwfYb648C3SngFGC5JBW2jCh5vB/5PRAwq+NUvIh7u5NyRRf0Xb/F6uP62Ay8VvTYwIi7soO524EDR9UZ1cFyuiIhVETETOIHcFFZLJ/V1VndHiq99cDroDaBfwWtDu9D3TuADkgYW9b2jhHqsjjjwrRJWkgvOL0rqI+lTwOSC1/8ZuE7SR5TTX9In8oG0EngXuEFSb0kzi87tyOH6ew7Yn7/Je4ykXpJOlXRmcSf56YzHgLmS+kkaB1zV0QUlHS3pCknHRsQ7wO/ITTUBvAYMlnRsaX9c7/HV/LXHA38NfD/fvg64UNIHJA0FvlR03mtAhz8fEBHbgf8C5klqyN/QvgYovmFsdc6Bb4mLiD8AnwI+D/wGmEUuSA++vhqYTe7G4G+BLfljC8+9BtgLXAn8iNycfGfXO1x/7wLTgYnAS8BucvPXnYXxDcAAclMci8ivAurE54Bt+VU31wFX5K/5PLmbzVvz00hdmZZ5Nl//T4B/jIhl+fbvkbvxvA1Yxp++EBw0D7gtf72beL/LgEZyo/3Hyd0D+XEX6rI6IH8AiqWdpJ8B90bE4cLXzI7AI3xLHUn/Q9LQ/JTOVcBfAP9W67rMejrfhbc0OoXcDdD+wFbgMxHxSm1LMuv5PKVjZpYRntIxM8uIVE/pDBkyJBobG2tdhplZj7JmzZrdEXF8cXsqA1/SDGDGmDFjWL16da3LMTPrUSR1+BPiqZzSiYgnImLOscd25+dWzMysI6kMfEkzJC3Yt29frUsxM6sbqQx8j/DNzJKX+jl8M7OueOedd2hra+Ott9468sE9XENDAyNGjKBPnz4lHZ/qdfhNTU3hm7Zm1hUvvfQSAwcOZPDgwbx309X6EhHs2bOH/fv3M3r06Pe8JmlNRLzvU8tSOaVjZtZdb731Vt2HPYAkBg8e3KXvZFIZ+L5pa2blqPewP6ir7zOVge+btmZmyUvlTVs7vMbmJ6t2rW23f6Jq1zKrhKT/v5Tyf2Lv3r0sXryY66+/vkt9X3jhhSxevJhBgwZ1t7zDSmXge5VO9wwc25x4nxPuP3Kf669an/h1zXqyvXv3cs8997wv8A8cOEDv3p3Hbmtra0Xr8pSOmVnCmpubefHFF5k4cSJnnnkmZ599NhdddBHjxo0D4OKLL+aMM85g/PjxLFiw4NB5jY2N7N69m23btjF27Fhmz57N+PHj+fjHP86bb75Zdl2pDHwzs57s9ttv50Mf+hDr1q3jW9/6FmvXruXb3/42v/rVrwBYuHAha9asYfXq1dx1113s2bPnfX1s3ryZL3zhC2zcuJFBgwbx6KOPll1XKgPfq3TMrJ5Mnjz5PWvl77rrLk477TSmTJnC9u3b2bx58/vOGT16NBMnTgTgjDPOYNu2bWXXkcrA95SOmdWT/v37H3r8zDPP8OMf/5iVK1fyi1/8gkmTJnW4lr5v376HHvfq1YsDBw6UXUcqA9/MrCcbOHAg+/fv7/C1ffv2cdxxx9GvXz+ef/55fvrTn1atrlSu0rGeJcllb14Gakmrxb+pwYMHM3XqVE499VSOOeYYTjzxxEOvTZs2jXvvvZexY8dyyimnMGXKlKrV5cA3M6uAxYsXd9jet29fnnrqqQ5fOzhPP2TIEDZs2HCo/aabbkqkplQGvtfh9yxJrv8vZd1/If8MgFnpUjmH75u2ZmbJS2Xgm5lZ8hz4ZmYZ4cA3M8uIVN607YmquYOlmVl3OPDNrL7NTXjxx9zkt3wZMGAAr7/+euL9FkvllI730jEzS14qA9/LMs2sJ2tubmb+/PmHns+dO5dvfvObnHvuuZx++ulMmDCBH/7wh1WvK5WBb2bWk82aNYuWlpZDz1taWrjqqqt4/PHHWbt2LcuXL+fGG28kIqpal+fwK6ASnzxlZj3HpEmT2LVrFzt37qS9vZ3jjjuOoUOH8uUvf5kVK1Zw1FFHsWPHDl577TWGDh1atboc+GZmFfDZz36WJUuW8OqrrzJr1iweeugh2tvbWbNmDX369KGxsbHDbZEryYFvZlYBs2bNYvbs2ezevZtnn32WlpYWTjjhBPr06cPy5ct5+eWXq16TA9/M6lsFllGWYvz48ezfv5/hw4dz0kknccUVVzBjxgwmTJhAU1MTH/7wh6tekwPfzKxC1q//026uQ4YMYeXKlR0eV401+OBVOmZmmVG1wJf055Luk7SkWtc0M7M/KSnwJS2UtEvShqL2aZJekLRF0mHXIkbE1oi4ppxizcys+0qdw18E3A08cLBBUi9gPnAe0AaskrQU6AXMKzr/6ojYVXa1ZmbWbSUFfkSskNRY1DwZ2BIRWwEkPQLMjIh5wPTuFiRpDjAHYNSoUd3txszMipQzhz8c2F7wvC3f1iFJgyXdC0ySdGtnx0XEgohoioim448/vozyzMysUNWWZUbEHuC6Uo71h5ibWVIm3D8h0f7WX7X+iMfs3buXxYsXc/3113e5/zvvvJM5c+bQr1+/7pR3WOWM8HcAIwuej8i3lc27ZZpZT7Z3717uueeebp1755138vvf/z7hinLKGeGvAk6WNJpc0F8KXJ5EUR7hm1lP1tzczIsvvsjEiRM577zzOOGEE2hpaeHtt9/mk5/8JN/4xjd44403uOSSS2hra+Pdd9/lq1/9Kq+99ho7d+7kox/9KEOGDGH58uWJ1lVS4Et6GDgHGCKpDfh6RNwn6QbgaXIrcxZGxMYkioqIJ4AnmpqaZifRn5lZNd1+++1s2LCBdevWsWzZMpYsWcJzzz1HRHDRRRexYsUK2tvbGTZsGE8+mft41H379nHsscdyxx13sHz5coYMGZJ4XaWu0rmsk/ZWoDXRivAI38zqx7Jly1i2bBmTJk0CctsobN68mbPPPpsbb7yRW265henTp3P22WdXvJZU7qXjEb6Z1YuI4NZbb+Xaa69932tr166ltbWV2267jXPPPZevfe1rFa3Fe+mYmSVs4MCB7N+/H4Dzzz+fhQsXHtogbceOHYc+HKVfv35ceeWV3Hzzzaxdu/Z95yYtlSN8T+mYWVJKWUaZtMGDBzN16lROPfVULrjgAi6//HLOOussAAYMGMCDDz7Ili1buPnmmznqqKPo06cP3/nOdwCYM2cO06ZNY9iwYYnftFW1P1OxK5qammL16tW1LqMkjc1PHnrsjzisnlr8Z7Z027RpE2PHjq11GVXT0fuVtCYimoqP9ZSOmVlGpDLwJc2QtGDfvtp8Uo2ZWT1KZeD7J23NrBxpnqpOUlffZyoD38ysuxoaGtizZ0/dh35EsGfPHhoaGko+x6t0zKyujBgxgra2Ntrb22tdSsU1NDQwYsSIko9PZeD7B6/MrLv69OnD6NGja11GKnlKx8wsIxz4ZmYZ4cA3M8uIVAa+1+GbmSUvlYHvdfhmZslL5Sods1Il/XmlpfIePtYTpXKEb2ZmyXPgm5llhAPfzCwjHPhmZhmRysD3skwzs+SlMvC9LNPMLHmpDHwzM0ueA9/MLCMc+GZmGeHANzPLCAe+mVlGOPDNzDKiqpunSboY+ATwZ8B9EbGsmtc3M8uykkf4khZK2iVpQ1H7NEkvSNoiqflwfUTEDyJiNnAdMKt7JZuZWXd0ZYS/CLgbeOBgg6RewHzgPKANWCVpKdALmFd0/tURsSv/+Lb8eWZmViUlB35ErJDUWNQ8GdgSEVsBJD0CzIyIecD04j4kCbgdeCoi1na3aDMz67pyb9oOB7YXPG/Lt3Xmb4GPAZ+RdF1HB0iaI2m1pNXt7e1llmdmZgdV9aZtRNwF3HWEYxYACwCampqiGnWZmWVBuSP8HcDIgucj8m1l8W6ZZmbJKzfwVwEnSxot6WjgUmBp+WWZmVnSurIs82FgJXCKpDZJ10TEAeAG4GlgE9ASERvLLcrbI5uZJa8rq3Qu66S9FWhNrCIzM6uIVG6t4Dl8M7PkpTLwPaVjZpa8VAa+R/hmZslLZeB7hG9mlrxUBr6ZmSUvlYHvKR0zs+SlMvA9pWNmlrxUBr6ZmSXPgW9mlhGpDHzP4ZuZJS+Vge85fDOz5KUy8M3MLHkOfDOzjHDgm5llRCoD3zdtzcySl8rA901bM7PkVfVDzKupsfnJWpdgZpYqqRzhm5lZ8up2hL+t4fKaXXsCo2p2bTOzzniEb2aWEQ58M7OMSGXge1mmmVnyUhn4XpZpZpa8VAa+mZklz4FvZpYRDnwzs4xw4JuZZYQD38wsIxz4ZmYZUbXAlzRW0r2Slkj6n9W6rpmZ5ZQU+JIWStolaUNR+zRJL0jaIqn5cH1ExKaIuA64BJja/ZLNzKw7Sh3hLwKmFTZI6gXMBy4AxgGXSRonaYKkHxX9OiF/zkXAk0BrYu/AzMxKUtJumRGxQlJjUfNkYEtEbAWQ9AgwMyLmAdM76WcpsFTSk8Di7hZtZmZdV872yMOB7QXP24CPdHawpHOATwF9OcwIX9IcYA7AqFHeZtjMLClV2w8/Ip4BninhuAWSXgFmHH300WdUui4zs6woJ/B3ACMLno/It5UtIp4AnmhqapqdRH9mdWVujTYVnOvda3u6cpZlrgJOljRa0tHApcDSJIry9shmZskrdVnmw8BK4BRJbZKuiYgDwA3A08AmoCUiNiZRlLdHNjNLXqmrdC7rpL0VL7E0M+sRUrm1gqd0zMySl8rA95SOmVnyUhn4HuGbmSUvlYHvEb6ZWfJSGfhmZpa8VAa+p3TMzJKXysD3lI6ZWfJSGfhmZpY8B76ZWUakMvA9h29mlrxUBr7n8M3MkpfKwDczs+Q58M3MMsKBb2aWEakMfN+0NTNLXioD3zdtzcySV7UPMTerK7X6XFlgwuhRNbnu+ppc1ZKUyhG+mZklz4FvZpYRDnwzs4xw4JuZZUQqA9/LMs3MkpfKwPeyTDOz5NXtssxaLV0zM0urVI7wzcwseQ58M7OMqNspHbNK8pSh9UQe4ZuZZYQD38wsI6oa+JL6S1otaXo1r2tmZiUGvqSFknZJ2lDUPk3SC5K2SGouoatbgJbuFGpmZuUp9abtIuBu4IGDDZJ6AfOB84A2YJWkpUAvYF7R+VcDpwG/BBrKK9nMzLqjpMCPiBWSGouaJwNbImIrgKRHgJkRMQ9435SNpHOA/sA44E1JrRHxxw6OmwPMARg1yishzKyGavW5B3Mrs61MOcsyhwPbC563AR/p7OCI+F8Akj4P7O4o7PPHLQAWADQ1NUUZ9ZmZWYGqr8OPiEVHOkbSDGDGmDFjKl+QmVlGlLNKZwcwsuD5iHxb2bx5mplZ8soZ4a8CTpY0mlzQXwpcnkRRHuGb2XvU8DOE60mpyzIfBlYCp0hqk3RNRBwAbgCeBjYBLRGxMYmiPMI3M0teqat0LuukvRVoTbQiPMI3SyWPsnu8VG6t4BG+mVnyUhn4ZmaWvFQGvj/T1swseakMfE/pmJklL5WBb2ZmyUtl4HtKx8wseakMfE/pmJklL5WBb2ZmyXPgm5llRCoD33P4ZmbJS2Xgew7fzCx5qQx8MzNLngPfzCwjHPhmZhmRysD3TVszs+SlMvB909bMLHmpDHwzM0ueA9/MLCMc+GZmGeHANzPLCAe+mVlG9K51AR2RNAOYMWbMmFqXYmYZNmH0qJpcd32F+k3lCN/LMs3MkpfKwDczs+Q58M3MMsKBb2aWEQ58M7OMcOCbmWWEA9/MLCOqFviSzpH0H5LulXROta5rZmY5JQW+pIWSdknaUNQ+TdILkrZIaj5CNwG8DjQAbd0r18zMuqvUn7RdBNwNPHCwQVIvYD5wHrkAXyVpKdALmFd0/tXAf0TEs5JOBO4AriivdDMz64qSAj8iVkhqLGqeDGyJiK0Akh4BZkbEPGD6Ybr7LdC3sxclzQHmAIwaVZsfazYzq0fl7KUzHNhe8LwN+EhnB0v6FHA+MIjcdwsdiogFwAKApqamKKM+M0tQrfaVseRUbfO0iHgMeKyUY715mplZ8spZpbMDGFnwfES+rWzePM3MLHnlBP4q4GRJoyUdDVwKLE2iKEkzJC3Yt29fEt2ZmRmlL8t8GFgJnCKpTdI1EXEAuAF4GtgEtETExiSK8gjfzCx5pa7SuayT9lagNdGK8By+mVklpHJrBY/wzcySl8rANzOz5KUy8H3T1swseakMfE/pmJklL5WBb2ZmyUtl4HtKx8wseakMfE/pmJklL5WBb2ZmyXPgm5llRCoD33P4ZmbJS2Xgew7fzCx5qQx8MzNLngPfzCwjHPhmZhmRysD3TVszs+SlMvB909bMLHmpDHwzM0ueA9/MLCMUEbWuoVOS2oGXu3n6EGB3guX0BH7P2eD3XP/Kfb8fjIjjixtTHfjlkLQ6IppqXUc1+T1ng99z/avU+/WUjplZRjjwzcwyop4Df0GtC6gBv+ds8HuufxV5v3U7h29mZu9VzyN8MzMr4MA3M8uIugx8SdMkvSBpi6TmWtdTaZJGSlou6ZeSNkr6u1rXVA2Sekn6uaQf1bqWapA0SNISSc9L2iTprFrXVGmSvpz/N71B0sOSGmpdU9IkLZS0S9KGgrYPSPp/kjbnfz8uiWvVXeBL6gXMBy4AxgGXSRpX26oq7gBwY0SMA6YAX8jAewb4O2BTrYuoom8D/xYRHwZOo87fu6ThwBeBpog4FegFXFrbqipiETCtqK0Z+ElEnAz8JP+8bHUX+MBkYEtEbI2IPwCPADNrXFNFRcQrEbE2/3g/uSAYXtuqKkvSCOATwHdrXUs1SDoW+O/AfQAR8YeI2FvbqqqiN3CMpN5AP2BnjetJXESsAH5T1DwTuD//+H7g4iSuVY+BPxzYXvC8jToPv0KSGoFJwM9qW0nF3Qn8PfDHWhdSJaOBduBf8tNY35XUv9ZFVVJE7AD+Efg18AqwLyKW1baqqjkxIl7JP34VODGJTusx8DNL0gDgUeBLEfG7WtdTKZKmA7siYk2ta6mi3sDpwHciYhLwBgl9m59W+XnrmeS+2A0D+ku6srZVVV/k1s4nsn6+HgN/BzCy4PmIfFtdk9SHXNg/FBGP1bqeCpsKXCRpG7kpu7+U9GBtS6q4NqAtIg5+57aE3BeAevYx4KWIaI+Id4DHgP9W45qq5TVJJwHkf9+VRKf1GPirgJMljZZ0NLmbPEtrXFNFSRK5ud1NEXFHreuptIi4NSJGREQjub/ff4+Iuh75RcSrwHZJp+SbzgV+WcOSquHXwBRJ/fL/xs+lzm9UF1gKXJV/fBXwwyQ67Z1EJ2kSEQck3QA8Te6u/sKI2FjjsiptKvA5YL2kdfm2r0REaw1rsuT9LfBQfiCzFfjrGtdTURHxM0lLgLXkVqL9nDrcYkHSw8A5wBBJbcDXgduBFknXkNsi/pJEruWtFczMsqEep3TMzKwDDnwzs4xw4JuZZYQD38wsIxz4ZmYZ4cA3M8sIB76ZWUb8fwznX7tp92LtAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEICAYAAABcVE8dAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAaEUlEQVR4nO3df5RU5Z3n8fdHaGwRIiytMfxKEzEGhBG0h+iakzHHmICKJJlRJLqrE4aOUWedrHHEHZMlO9mVHWdcZdUYJnI0G1F7SGIgtkeSHQzZESM/RgMEVGQ1NKggIyzqYMR894+6kKLtbqq7blfd6vt5ncOh6tat5/ne7ubTTz334V5FBGZm1vcdVe0CzMysMhz4ZmY54cA3M8sJB76ZWU448M3McsKBb2aWEw58qyhJV0r6P9WuoxSSnpD0Z8njyyQtT7HtjZLOSR7Pk/T9FNv+T5K+m1Z71nc48K0mSQpJYyvVX0Q8EBGfOdJ+ku6T9K0S2js1Ip4oty5J50hqa9f2f4uIPyu3bet7HPiWO5L657FvMwe+9QpJoyT9UNIuSbsl3dnBPo3JSL1/0bbiaZSxkn4uaa+k1yU9nGxfmez+rKQ3Jc1Mtl8o6RlJeyQ9KekPitp9SdKNkn4FvNVR8Eo6T9LmpL87ARW9dmgqSgX/Q9JOSf9P0npJEyQ1A5cBf5nUtayzvpNtny7qvl7Sw5L2SVon6bSivg/7NHPwU4SkY4HHgOFJf29KGt5+ikjSRckU0p7k6zuu3dfla5J+lRz3w5Lqu/zmWs1y4FvqJPUDfgK8DDQCI4CHetDUXwPLgaHASOB/AkTEJ5PXT4uIQRHxsKTJwCLgy8Aw4DvAUklHF7U3C7gAGBIRB9rV3AD8ELgZaABeBM7upK7PAJ8EPgocB1wC7I6IhcADwN8kdU0vpe/EDOAfgH8DLAYekVTX1RcnIt4CpgE7kv4GRcSOdsf1UeBB4C+A44FWYJmkAUW7XQJMBcYAfwBc2VW/Vrsc+NYbpgDDgRsi4q2I2B8RPTlR+y7wYWB4CW00A9+JiF9GxHsRcT/wDnBm0T4LImJbRPxrB+8/H9gYEUsi4l3gduDVLuoaDHwMUERsiohXjnAsXfUNsLao79uA+na199RM4NGI+GnS9t8CxwD/tl1tOyLiX4BlwKQU+rUMcuBbbxgFvNzJSLY7/pLCtMrTyZTEl7rY98PA9cm0xR5Je5I6hhfts62L9w8vfj0KVxXscP+I+EfgTuAuYKekhZI+cIRj6arvw16PiN8BbRxee08Np/BJq7jtbRQ+dR1U/IvtbWBQCv1aBjnwrTdsA0aXcILyreTvgUXbTjz4ICJejYg5ETGcwlTN3V2szNkG/NeIGFL0Z2BEPFi0T1eXhn2Fwi8IoDBPX/y8vYhYEBFnAOMpTO3ccIQ+jnRZ2uK+j6IwhXVweuZtOvkaldDuDgq/DA+2ffC4th/hfdYHOfCtNzxNIUDnSzpWUr2k982HR8QuCsFzuaR+yQj+pIOvS7pY0sjk6RsUwu13yfPXgI8UNff3wFWSPp6cVD1W0gWSBpdY86PAqZK+kPyi+g8cHqyHSPrDpJ86Cr+09ndRV6nOKOr7LyhMRz2VvPYM8MXkazQV+KOi970GDJN0XCfttgAXSDo3qff6pO0ne1Cj1TgHvqUuIt4DpgNjgd9QmJ6Y2cnucyiMjncDp3J4EP0h8EtJbwJLgesiYmvy2jzg/mT65pKIWJO0dSeFXw5b6MbJx4h4HbgYmJ/UcjLwT53s/gEKv2DeoDBdshu4NXntXmB8UtcjpfYP/JjC1+gN4N8BX0jm3AGuo/D13ENhFdChdiNiM4WTsluTPg+bBoqI54DLKZzwfj1pZ3pE/LYbtVkfId8AxcwsHzzCNzPLCQe+mVlOOPDNzHLCgW9mlhOZvpBTQ0NDNDY2VrsMM7Oasnbt2tcj4vj22zMd+I2NjaxZs6baZZiZ1RRJL3e03VM6ZmY54cA3M8sJB76ZWU5keg7fzKy73n33Xdra2ti/f3+1S+l19fX1jBw5krq6Lm+dcIgD38z6lLa2NgYPHkxjYyOFi4P2TRHB7t27aWtrY8yYMSW9p2JTOpI+IuleSUsq1aeZ5c/+/fsZNmxYnw57AEkMGzasW59kygp8SYuS+3puaLd9qqTnJG2RNBcgIrZGxOxy+jMzK0VfD/uDunuc5Y7w76NwL8ziAvpRuBPQNAo3h5glaXyZ/ZiZWZnKmsOPiJWSGtttngJsOXjdckkPUbhB869LaVNSM4X7kzJ69Oge1zbx/ok9fm+51l+xvmp9m9nhGuc+mmp7L82/4Ij77Nmzh8WLF3P11Vd3q+3zzz+fxYsXM2TIkJ6W16XemMMfweH372wDRkgaJukeYLKkmzp7c0QsjIimiGg6/vj3/c9gM7PM27NnD3fffff7th840PVtnltbW3st7KGCq3QiYjdwVSn7SpoOTB87trPbl5qZZdfcuXN58cUXmTRpEnV1ddTX1zN06FA2b97M888/z+c+9zm2bdvG/v37ue6662hubgZ+fzmZN998k2nTpvGJT3yCJ598khEjRvDjH/+YY445pqy6emOEv53Db/48km7eMDkilkVE83HHdXabTjOz7Jo/fz4nnXQSzzzzDLfeeivr1q3jjjvu4Pnnnwdg0aJFrF27ljVr1rBgwQJ27979vjZeeOEFrrnmGjZu3MiQIUP4wQ9+UHZdvRH4q4GTJY2RNAC4lML9SEsmabqkhXv37u2F8szMKmvKlCmHrZVfsGABp512GmeeeSbbtm3jhRdeeN97xowZw6RJkwA444wzeOmll8quo9xlmQ8Cq4BTJLVJmh0RB4BrgceBTUBLRGzsTrse4ZtZX3LsscceevzEE0/ws5/9jFWrVvHss88yefLkDtfSH3300Yce9+vX74jz/6Uod5XOrE62twKtPW3Xc/hmVssGDx7Mvn37Onxt7969DB06lIEDB7J582aeeuqpitWVyUsrRMQyYFlTU9OcatdiZrWtlGWUaRs2bBhnn302EyZM4JhjjuGDH/zgodemTp3KPffcw7hx4zjllFM488wzK1ZXJgPfI3wzq3WLFy/ucPvRRx/NY4891uFrB+fpGxoa2LDh9xcw+NrXvpZKTZm8PLLn8M3M0pfJwDczs/RlMvC9LNPMLH2ZDHxP6ZiZpS+TgW9mZunLZOB7SsfMLH2ZXJbpdfhmlpp5KU8Nz0t/IDpo0CDefPPN1NttL5MjfDMzS18mR/hmZrVs7ty5jBo1imuuuQaAefPm0b9/f1asWMEbb7zBu+++y7e+9S1mzJhR0bo8wjczS9nMmTNpaWk59LylpYUrrriCH/3oR6xbt44VK1Zw/fXXExEVrSuTI3xfWsHMatnkyZPZuXMnO3bsYNeuXQwdOpQTTzyRr371q6xcuZKjjjqK7du389prr3HiiSdWrK5MBr5P2ppZrbv44otZsmQJr776KjNnzuSBBx5g165drF27lrq6OhobGzu8LHJvymTgm5nVupkzZzJnzhxef/11fv7zn9PS0sIJJ5xAXV0dK1as4OWXX654TQ58M+vbemEZZSlOPfVU9u3bx4gRI/jQhz7EZZddxvTp05k4cSJNTU187GMfq3hNDnwzs16yfv36Q48bGhpYtWpVh/tVYg0+eJWOmVluZDLwfWkFM7P0ZTLwfbVMM7P0ZTLwzcwsfQ58M7OccOCbmeWEl2WaWZ828f6Jqba3/or1R9xnz549LF68mKuvvrrb7d9+++00NzczcODAnpTXJY/wzcxStmfPHu6+++4evff222/n7bffTrmigoqN8CUdC9wN/BZ4IiIeqFTfZmaVNHfuXF588UUmTZrEeeedxwknnEBLSwvvvPMOn//85/nmN7/JW2+9xSWXXEJbWxvvvfceX//613nttdfYsWMHn/rUp2hoaGDFihWp1lVW4EtaBFwI7IyICUXbpwJ3AP2A70bEfOALwJKIWCbpYcCBb2Z90vz589mwYQPPPPMMy5cvZ8mSJTz99NNEBBdddBErV65k165dDB8+nEcffRSAvXv3ctxxx3HbbbexYsUKGhoaUq+r3Cmd+4CpxRsk9QPuAqYB44FZksYDI4FtyW7vldmvmVlNWL58OcuXL2fy5MmcfvrpbN68mRdeeIGJEyfy05/+lBtvvJFf/OIXVOL/HZU1wo+IlZIa222eAmyJiK0Akh4CZgBtFEL/GXzuwMxyIiK46aab+PKXv/y+19atW0drays333wz5557Lt/4xjd6tZbeCN4R/H4kD4WgHwH8EPhjSd8GlnX2ZknNktZIWrNr165eKM/MrHcNHjyYffv2AfDZz36WRYsWHbpA2vbt2w/dHGXgwIFcfvnl3HDDDaxbt+59701bxU7aRsRbwJ+WsN9CSa8A0wcMGHBG71dmZn1ZKcso0zZs2DDOPvtsJkyYwLRp0/jiF7/IWWedBcCgQYP4/ve/z5YtW7jhhhs46qijqKur49vf/jYAzc3NTJ06leHDh6d+0lbl3lMxmdL5ycGTtpLOAuZFxGeT5zcBRMQt3W27qakp1qxZ06O60l572x3V+AEzs4JNmzYxbty4apdRMR0dr6S1EdHUft/emNJZDZwsaYykAcClwNLuNOCrZZqZpa+swJf0ILAKOEVSm6TZEXEAuBZ4HNgEtETExu6066tlmpmlr9xVOrM62d4KtPa0XUnTgeljx47taRNmlmMRgaRql9Hrujsln8nlkR7hm1lP1dfXs3v37m6HYa2JCHbv3k19fX3J78nkxdM8wjeznho5ciRtbW3kYVl3fX09I0eOLHn/TAZ+RCwDljU1Nc2pdi1mVlvq6uoYM2ZMtcvIpExO6ZiZWfoyGfhelmlmlr5MBr5P2pqZpS+TgW9mZulz4JuZ5UQmA99z+GZm6ctk4HsO38wsfZkMfDMzS58D38wsJzIZ+J7DNzNLXyYD33P4Zmbpy2Tgm5lZ+hz4ZmY54cA3M8sJB76ZWU5kMvC9SsfMLH2ZDHyv0jEzS18mA9/MzNLnwDczywkHvplZTjjwzcxywoFvZpYTDnwzs5yoWOBL+oikeyUtqVSfZmb2eyUFvqRFknZK2tBu+1RJz0naImluV21ExNaImF1OsWZm1nP9S9zvPuBO4HsHN0jqB9wFnAe0AaslLQX6Abe0e/+XImJn2dWamVmPlRT4EbFSUmO7zVOALRGxFUDSQ8CMiLgFuLCnBUlqBpoBRo8e3dNmzMysnXLm8EcA24qetyXbOiRpmKR7gMmSbupsv4hYCHwTWDdgwIAyyjMzs2IVO2kbEbsj4qqIOCn5FNDVvr6WjplZysoJ/O3AqKLnI5NtZfPVMs3M0ldO4K8GTpY0RtIA4FJgaRpFeYRvZpa+UpdlPgisAk6R1CZpdkQcAK4FHgc2AS0RsTGNojzCNzNLX6mrdGZ1sr0VaE21okK7y4BlTU1Nc9Ju28wsr3xpBTOznMhk4HtKx8wsfZkMfJ+0NTNLXyYD3yN8M7P0ZTLwPcI3M0tfJgPfzMzS58A3M8uJTAa+5/DNzNKXycD3HL6ZWfoyGfhmZpY+B76ZWU5kMvA9h29mlr5S72lbUb54Wm2ZeP/EqvW9/or1Vevb+r5q/Wz31s91Jkf4ZmaWPge+mVlOOPDNzHLCgW9mlhMOfDOznMhk4HtZpplZ+jIZ+L60gplZ+jIZ+GZmlj4HvplZTjjwzcxywoFvZpYTDnwzs5yo6MXTJH0OuAD4AHBvRCyvZP9mZnlWcuBLWgRcCOyMiAlF26cCdwD9gO9GxPzO2oiIR4BHJA0F/hbok4Hf166wZ2Z9Q3dG+PcBdwLfO7hBUj/gLuA8oA1YLWkphfC/pd37vxQRO5PHNyfvMzOzCik58CNipaTGdpunAFsiYiuApIeAGRFxC4VPA4eRJGA+8FhErOuoH0nNQDPA6NGjSy3PzMyOoNyTtiOAbUXP25Jtnflz4NPAn0i6qqMdImJhRDRFRNPxxx9fZnlmZnZQRU/aRsQCYMGR9pM0HZg+duzY3i/KrNbMq9IlR+b52la1rtwR/nZgVNHzkcm2svhaOmZm6Ss38FcDJ0saI2kAcCmwtNyifLVMM7P0lRz4kh4EVgGnSGqTNDsiDgDXAo8Dm4CWiNhYblEe4ZuZpa87q3RmdbK9FWhNrSI8h29m1hsyeWkFj/DNzNKXycD3HL6ZWfoyGfge4ZuZpS+TgW9mZunLZOB7SsfMLH2ZDHxP6ZiZpS+TgW9mZunLZOB7SsfMLH2ZDHxP6ZiZpS+TgW9mZulz4JuZ5YQD38wsJzIZ+D5pa2aWvkwGvk/ampmlL5OBb2Zm6XPgm5nlREVvYm5m5Zs4ZnRV+l1flV4tTR7hm5nlRCYD36t0zMzSl8nA9yodM7P0ZTLwzcwsfQ58M7OccOCbmeWEA9/MLCcc+GZmOVGxwJc0TtI9kpZI+kql+jUzs4KSAl/SIkk7JW1ot32qpOckbZE0t6s2ImJTRFwFXAKc3fOSzcysJ0od4d8HTC3eIKkfcBcwDRgPzJI0XtJEST9p9+eE5D0XAY8CrakdgZmZlaSka+lExEpJje02TwG2RMRWAEkPATMi4hbgwk7aWQoslfQosLijfSQ1A80Ao0dX55ohZmZ9UTkXTxsBbCt63gZ8vLOdJZ0DfAE4mi5G+BGxEFgI0NTUFGXUZ2ZmRSp2tcyIeAJ4opR9JU0Hpo8dO7Y3SzIzy5VyVulsB0YVPR+ZbCubr6VjZpa+cgJ/NXCypDGSBgCXAkvTKMpXyzQzS1+pyzIfBFYBp0hqkzQ7Ig4A1wKPA5uAlojYmEZRHuGbmaWv1FU6szrZ3kovLLH0HL6ZWfoyeWkFj/DNzNKXycD3HL6ZWfoyGfge4ZuZpS+TgW9mZunLZOB7SsfMLH0V+5+23RERy4BlTU1Nc6pdSy2ZeP/EapdgZhmWyRG+mZmlz4FvZpYTmQx8z+GbmaUvk4HvZZlmZunLZOCbmVn6HPhmZjmRycD3HL6ZWfoyGfiewzczS18mA9/MzNLnwDczywkHvplZTjjwzcxyIpOB71U6Zmbpy2Tge5WOmVn6Mhn4ZmaWPge+mVlOOPDNzHLCgW9mlhMOfDOznHDgm5nlREUDX9KxktZIurCS/ZqZWYmBL2mRpJ2SNrTbPlXSc5K2SJpbQlM3Ai09KdTMzMrTv8T97gPuBL53cIOkfsBdwHlAG7Ba0lKgH3BLu/d/CTgN+DVQX17JZmbWEyUFfkSslNTYbvMUYEtEbAWQ9BAwIyJuAd43ZSPpHOBYYDzwr5JaI+J3HezXDDQDjB49uuQDMTOzrpU6wu/ICGBb0fM24OOd7RwRfwUg6Urg9Y7CPtlvoaRXgOkDBgw4o4z6zMysSMVX6UTEfRHxkyPs42vpmJmlrJzA3w6MKno+MtlWNl8t08wsfeUE/mrgZEljJA0ALgWWplGUR/hmZukrdVnmg8Aq4BRJbZJmR8QB4FrgcWAT0BIRG9MoyiN8M7P0lbpKZ1Yn21uB1lQrKrS7DFjW1NQ0J+22zczyypdWMDPLiUwGvqd0zMzSl8nA90lbM7P0ZTLwPcI3M0tfJgPfI3wzs/RlMvDNzCx9Dnwzs5zIZOB7Dt/MLH2ZDHzP4ZuZpS+TgW9mZulz4JuZ5UQmA99z+GZm6ctk4HsO38wsfZkMfDMzS18597Q1y63GuY9Wre/B46rWtdU4j/DNzHLCgW9mlhOZDHyv0jEzS18mA9+rdMzM0pfJwDczs/Q58M3McsKBb2aWEw58M7OccOCbmeWEIqLaNXRK0i7g5R6+vQF4PcVyaoGPOR98zH1fucf74Yg4vv3GTAd+OSStiYimatdRST7mfPAx9329dbye0jEzywkHvplZTvTlwF9Y7QKqwMecDz7mvq9XjrfPzuGbmdnh+vII38zMijjwzcxyouYDX9JUSc9J2iJpbgevHy3p4eT1X0pqrHyV6SrhmP+jpF9L+pWk/y3pw9WoM01HOuai/f5YUkiq6SV8pRyvpEuS7/NGSYsrXWPaSvi5Hi1phaR/Tn62z69GnWmStEjSTkkbOnldkhYkX5NfSTq9rA4jomb/AP2AF4GPAAOAZ4Hx7fa5GrgneXwp8HC1667AMX8KGJg8/koejjnZbzCwEngKaKp23b38PT4Z+GdgaPL8hGrXXYFjXgh8JXk8Hnip2nWncNyfBE4HNnTy+vnAY4CAM4FfltNfrY/wpwBbImJrRPwWeAiY0W6fGcD9yeMlwLmSVMEa03bEY46IFRHxdvL0KWBkhWtMWynfZ4C/Bv47sL+SxfWCUo53DnBXRLwBEBE7K1xj2ko55gA+kDw+DthRwfp6RUSsBP6li11mAN+LgqeAIZI+1NP+aj3wRwDbip63Jds63CciDgB7gWEVqa53lHLMxWZTGCHUsiMec/JRd1REVO/u4ukp5Xv8UeCjkv5J0lOSplasut5RyjHPAy6X1Aa0An9emdKqqrv/3rvUv+xyLLMkXQ40AX9U7Vp6k6SjgNuAK6tcSiX1pzCtcw6FT3ArJU2MiD1Vrap3zQLui4i/k3QW8L8kTYiI31W7sFpR6yP87cCooucjk20d7iOpP4WPgrsrUl3vKOWYkfRp4K+AiyLinQrV1luOdMyDgQnAE5JeojDXubSGT9yW8j1uA5ZGxLsR8X+B5yn8AqhVpRzzbKAFICJWAfUULjLWl5X0771UtR74q4GTJY2RNIDCSdml7fZZClyRPP4T4B8jORtSo454zJImA9+hEPa1PrcLRzjmiNgbEQ0R0RgRjRTOW1wUEWuqU27ZSvm5foTC6B5JDRSmeLZWssiUlXLMvwHOBZA0jkLg76polZW3FPj3yWqdM4G9EfFKTxur6SmdiDgg6VrgcQpn+RdFxEZJ/wVYExFLgXspfPTbQuHkyKXVq7h8JR7zrcAg4B+S89O/iYiLqlZ0mUo85j6jxON9HPiMpF8D7wE3RETNfnIt8ZivB/5e0lcpnMC9ssYHb0h6kMIv7obk3MR/BuoAIuIeCucqzge2AG8Df1pWfzX+9TIzsxLV+pSOmZmVyIFvZpYTDnwzs5xw4JuZ5YQD38wsJxz4ZmY54cA3M8uJ/w+X6ujyZSncLgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_histogram('nodes', bins=50)\n",
    "plot_histogram('edges', bins=50)\n",
    "plot_histogram('degree')\n",
    "plot_histogram('cluster')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dbfe20cb",
   "metadata": {},
   "source": [
    "While the average number of nodes and edges are similar across the splits, the histograms show that the tails are slightly heavier for the training graphs."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1af88b23",
   "metadata": {},
   "source": [
    "## Plot molecules\n",
    "\n",
    "We borrow code from the PyTorch Geometric [GNN explanation example](https://colab.research.google.com/drive/1fLJbFPz0yMCQg81DdCP5I8jXw9LoggKO?usp=sharing#scrollTo=9Hh3YNASuYxm) to visualize molecules from the graph objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "24e61eea",
   "metadata": {},
   "outputs": [],
   "source": [
    "def draw_molecule(g, edge_mask=None, draw_edge_labels=False):\n",
    "    g = g.copy().to_undirected()\n",
    "    node_labels = {}\n",
    "    for u, data in g.nodes(data=True):\n",
    "        node_labels[u] = data['name']\n",
    "    pos = nx.planar_layout(g)\n",
    "    pos = nx.spring_layout(g, pos=pos)\n",
    "    if edge_mask is None:\n",
    "        edge_color = 'black'\n",
    "        widths = None\n",
    "    else:\n",
    "        edge_color = [edge_mask[(u, v)] for u, v in g.edges()]\n",
    "        widths = [x * 10 for x in edge_color]\n",
    "    nx.draw(g, pos=pos, labels=node_labels, width=widths,\n",
    "            edge_color=edge_color, edge_cmap=plt.cm.Blues,\n",
    "            node_color='azure')\n",
    "    \n",
    "    if draw_edge_labels and edge_mask is not None:\n",
    "        edge_labels = {k: ('%.2f' % v) for k, v in edge_mask.items()}    \n",
    "        nx.draw_networkx_edge_labels(g, pos, edge_labels=edge_labels,\n",
    "                                    font_color='red')\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "def to_molecule(data):\n",
    "    ATOM_MAP = ['H', 'He', 'Li', 'Be', 'B', 'C', 'N', 'O', 'F', 'Ne', 'Na', 'Mg', 'Al', 'Si', 'P',\n",
    "                'S', 'Cl', 'Ar', 'K', 'Ca', 'Sc', 'Ti', 'V', 'Cr', 'Mn', 'Fe', 'Co', 'Ni', 'Cu', 'Zn',\n",
    "                'Ga', 'Ge', 'As', 'Se', 'Br', 'Kr', 'Rb', 'Sr', 'Y']\n",
    "    g = to_networkx(data, node_attrs=['x'])\n",
    "    for u, data in g.nodes(data=True):\n",
    "        data['name'] = ATOM_MAP[data['x'][0]]\n",
    "        del data['x']\n",
    "    return g"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "00c10c7c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAt4AAAF2CAYAAABZM59BAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3de1yUZf7/8TcyMIAgKZQHtLbd2srW1EQN1DRdLfOQGNA3W8tsO+lWVn7N5Je11pa1ibll1u6WWqnfFE+loqupmwqmqWsnW1o7WJ4STBSGAYaZ3x8DbMAMzHCY4+v5ePBAb+778jMqzHuu+dzXFWKz2WwCAAAA0KJaebsAAAAAIBgQvAEAAAAPIHgDAAAAHkDwBgAAADyA4A0AAAB4AMEbAAAA8ACCNwAAAOABBG8AAADAAwjeAAAAgAcQvAEAAAAPIHgDAAAAHkDwBgAAADyA4A0AAAB4AMEbAAAA8ACCNwAAAOABBG8AAADAAwjeAAAAgAcYvF0AAADBzCqprPKzTVKI7LNi4WJ2DAg0ITabzebtIgAACDYWSaWVn50xSDKKWTIgUBC8AQDwsFJJZjfOj5A9gAPwb7yLBQCAB7kbulV5fmkL1ALAswjeAAB4iEXOQ/eKpUs1KDFRCdHRuqxjR6UOH67cnTurv25W/W0pAHwfwRsAAA9xNmv9SmamHp8yRY/MmKG8kyf12ZEjumvSJG1Yu9al6wH4B3q8AQDwAKukcw6OFxYWqmtCguYvXKgxaWkNjhMjZs0Af8X3LgAAHlDm5Pje3FyZzWaNTElp0jgAfB/BGwAAD7A6OX66oEBx8fEyGFxbNNDZOAB8H8EbAAAPcNbX2S4uTgX5+bJYXLt1kv5QwH8RvAEA8IAQJ8d7JyXJaDRq/Zo1TRoHgO8jeAMA4AHOnnBjY2P1+KxZmjp5statWSOTyaTy8nJtzs7WzGnTXB4HgO9jVRMAADzA2aomVZYvWaJX585V3qFDio6JUY9evfRoRob6JifXOI9VTQD/RfAGAMBDitW0TXAMklo3Uy0API8XzQAAeIjRy9cD8C6CNwAAHmKQFNHIayMqrwfgvwjeAAB4kFH/Dd9Wq2urckeI2W4gEBC8AQDwMKPsvdof7dypUrNZcnK7VVVPN6EbCAwEbwAAvOCbr77S8IED9c5rr8kYEqIw2YN2mOxBO0b20E17CRA4+H4GAMAL/vrXv8pgMOjWW25pdN83AP/CcoIAAHhYaWmpOnfurGuvvVYrV670djkAPIRWEwAAPGz16tXKz8/XPffc4+1SAHgQM94AAHjY4MGD9c033+jw4cNq1Yo5MCBY8N0OAIAH5eXladu2bbr77rsJ3UCQ4TseAAAP+tvf/iaDwaCJEyd6uxQAHkarCQAAHsJNlUBwY8YbAAAPqbqp8t577/V2KQC8gBlvAAA85LrrrtN3332n//znP/R3A0GI73oAADwgLy9P27dv56ZKIIjxnQ8AgAdU7VR55513ersUAF5CqwkAAC2stLRUCQkJGjRokLKysrxdDgAvYcYbAIAWtmrVKhUUFLBTJRDkmPEGAKCFcVMlAIkZbwAAmo1VklmSSVJx5efvf/xRn372GTdVAmDGGwCAprJIKq38XFt5WZkqKipkCAlRbESEDB6uDYDvIHgDANAEpbLPcrsqQpKxhWoB4Nt4zwsAgEZyN3Sr8vzSFqgFgO8jeAMA0AgWOQ/dK5Yu1aDERCVER+uyjh2VOny4cnfurP66WY7bUgAENoI3AACN4GzW+pXMTD0+ZYoemTFDeSdP6rMjR3TXpEnasHatS9cDCFz0eAMA4CarpHMOjhcWFqprQoLmL1yoMWlpDY4TI2bAgGDC9zsAAG4qc3J8b26uzGazRqakNGkcAIGJ4A0AgJusTo6fLihQXHy8DAbXFg10Ng6AwETwBgDATc56NNvFxakgP18Wi2u3TtLrCQQXgjcAAG4KcXK8d1KSjEaj1q9Z06RxAAQmgjcAAG5y9uQZGxurx2fN0tTJk7VuzRqZTCaVl5drc3a2Zk6b5vI4AAITq5oAAOAmZ6uaVFm+ZIlenTtXeYcOKTomRj169dKjGRnqm5xc4zxWNQGCC8EbAIBGKFbTNsExSGrdTLUA8A+80AYAoBGMXr4egP8heAMA0AgGSRGNvDai8noAwYXgDQBAIxn1s/DtYudmhJjtBoIVwRsAgCYwyt6rfSY/X+aSElU4WcO7qqeb0A0EL4I3AABNZJCU+cwzuvrSS2WoqFBY5bEw2YN2jOyhm/YSILjxMwAAgCay2WxavXq1Env1Uhsjc9oAHGPGGwCAJtq3b5++//57jR071tulAPBhBG8AAJpo1apVCg0N1ahRo7xdCgAfxgY6AAA0gc1m0+WXX64LL7xQmzdv9nY5AHwYM94AADTBoUOHlJeXp5SUFG+XAsDHEbwBAGiC1atXS5LGjBnj5UoA+DpaTQAAaIJevXrJaDQqJyfH26UA8HHMeAMA0Ejffvut9u/fz2omAFxC8AYAoJGq2kzo7wbgClpNAABopAEDBujs2bM6ePCgt0sB4AeY8QYAoBFOnjypXbt20WYCwGUEbwAAGmHt2rWy2WwEbwAuo9UEAIBGuOGGG3T48GHl5eUpJCTE2+UA8APMeAMA4KYzZ87ogw8+0NixYwndAFxG8AYAwE3r16+XxWJhNRMAbiF4AwDgplWrVqlTp07q06ePt0sB4EcI3gAAuMFkMik7O1spKSlq1YqnUQCu4ycGAABu2LRpk0pKSljNBIDbCN4AALhh1apVateuna699lpvlwLAzxC8AQBwUVlZmdatW6fRo0fLYDB4uxwAfobgDQCAi7Zv364zZ87QZgKgUQjeAAC4aNWqVWrdurWGDh3q7VIA+CGCNwAALqioqNCaNWs0YsQIRUREeLscAH6I4A0AgAtyc3N18uRJ2kwANBp3hgAA4IBVUlnlZ5skW2SkHp4+XcNHjPBuYQD8VojNZrN5uwgAAHyFRVJp5efaSktLZTQaZZBkFLNXANxD8AYAoFKpJLMb50fIHsABwBX0eAMAIPdDtyrPL22BWgAEJoI3ACDoWeQ4dK9YulSDEhOVEB2tyzp2VOrw4crdubPGOWY5bksBgNoI3gCAoOdo1vqVzEw9PmWKHpkxQ3knT+qzI0d016RJ2rB2rUvXA0Bt9HgDAIKaVdK5WscKCwvVNSFB8xcu1Ji0NJfGiRGzWQDqx88IAEBQK3NwbG9ursxms0ampDRpHAD4OYI3ACCoWR0cO11QoLj4eBkMri8Y6GgcAPg5gjcAIKg56rdsFxengvx8WSyu3zZJ3yaAhhC8AQBBLcTBsd5JSTIajVq/Zk2TxgGAnyN4AwCCmqMnwtjYWD0+a5amTp6sdWvWyGQyqby8XJuzszVz2rQ655vNZu368EN9//33LV8wAL/FqiYAgKBmlXTOZpNC6s5ZL1+yRK/Onau8Q4cUHROjHr166dGMDPVNTq5xXmlpqbp27qyC/HwlJycrPT1daWlp6tSpk4ceBQB/QPAGAAS1AwcO6Lv8fF07eLBCQ0MbNYZB0rGvvtLy5cu1fPlyffLJJwoJCVH//v2Vnp6u1NRUdejQoXkLB+B3CN4AgKBktVo1b948TZ8+XYOGDNGy999vdPBuLXv4rvLll19qxYoVevfdd/X5558rJCREAwcOVHp6um6++WZdcMEFzfIYAPgXgjcAIOicOHFCEyZM0KZNm3TTTTfpjTfeUHRcnMNt4xsSIclYz9c///zz6hD+5ZdfqlWrVrruuuuUnp6usWPHKj4+vpGPAoC/IXgDAILKhg0bdOedd+rs2bOaO3eu7r33XoVU9neXSm6F74ZC98/ZbDZ99tlnWr58ud5991199dVXCg0N1ZAhQ5Senq6UlBS1a9fOzUcDwJ8QvAEAQcFsNmv69OmaN2+eunXrpv/7v/9T165d65xnkT2A17eCt0H2wO369jo12Ww2HTx4sDqEf/311zIYDBo6dKjS09N10003qW3bto0cHYCvIngDAALeF198oXHjxungwYN68MEH9fzzzysiIqLea6yybwNvlX1znBDZlx4MV/OuxWuz2bR///7qGzO//fZbhYWF6frrr1d6erpGjx6t2NjYZvwTAXgLwRsAELBsNpv++te/6uGHH1Z0dLQWLlyoESNGeLssp2w2m/bu3Vsdwr///nuFh4frhhtu0C233KJRo0YpJibG22UCaCSCNwAgIBUUFOjuu+/W6tWrNWzYMC1atEgdO3b0dlkus1qt2rNnj959912tWLFCR48eldFo1I033qj09HSNHDlS0dHRzfNnyTOz+0CwI3gDAALOtm3bNH78eP3444+aPXu2pkyZolat/DdCWq1W5ebmavny5VqxYoWOHz+uyMhIjRgxQunp6RoxYoSioqLcHtcT/ewA/ovgDQAIGOXl5XryySc1e/Zs/frXv9ayZcvUs2dPb5fVrCoqKrRr1y4tX75cWVlZOnnypKKiojRq1Cilp6dr+PDhioyMbHCcllzBBYBjBG8AQEA4fPiwxo0bpz179uj3v/+9XnrpJbVu3drbZbWoiooKffjhh1q+fLlWrlypU6dOKTo6WqNHj1Z6erquv/56hzeRuhu6qxC+gaYheAMA/N7bb7+tSZMmyWAw6G9/+5tSU1O9XZLHWSwWbd++vTqEnz59Wm3atNFNN92k9PR0DR06VEajURZJxU7GWLF0qeZnZuqrL79UdEyMuvXooUczMpTUv3/1ObV36QTgOoI3AMBvFRYWavLkyVqyZIkGDBigd955RxdeeKG3y/K68vJybd26VcuXL9fq1av1008/KTY2VmPGjNHM2bMV1769VLlpUJVXMjP10uzZynztNQ25/nqFh4dry8aNyvnwQz395z9Xn2eQPXwDcB/BGwDgl3bv3q1x48bpyJEjevLJJzVjxgyFhoZ6uyyfU1ZWpi1btmj58uXasXOncj79VBG1esALCwvVNSFB8xcu1Ji0tAbHjBGrnQCNQfAGAPiViooKPffcc3rqqafUpUsXLVmyRMnJyd4uyy8UWSwqkxRqqNkssmXjRt0ycqROms0yGBpuJDHK3u8NwD20aQEA/Mb333+v3/3ud/rwww916623asGCBezq6IZWBoMcvSdwuqBAcfHxLoVuyb7eNwD38U4RAMAvrFy5Ut27d9f+/fu1ePFiLVmyhNDtJmdvcbeLi1NBfr4slvpW9G54HAD1I3gDAHxacXGx7rnnHqWmpuqSSy7RgQMHdPvttyuk1s2BaJizv7HeSUkyGo1av2ZNk8YBUD+CNwDAZx04cEC9evXS3//+d02fPl07d+7UJZdc4u2y/JazJ/3Y2Fg9PmuWpk6erHVr1shkMqm8vFybs7M1c9o0l8cBUD9urgQA+Byr1ap58+Zp+vTpio+P19tvv63Bgwd7uyy/Z5V0rp6vL1+yRK/Onau8Q4cUHROjHr166dGMDPWtdfMqq5oAjUPwBgD4lBMnTmjChAnatGmTbrrpJr3xxhuKi4vzdlkBo1iSa53cjrGON9B4vGAFAPiMDRs2qHv37vrnP/+pBQsWaPXq1YTuZtbULd/ZMh5oPII3AMDrzGazpkyZohEjRqh9+/bat2+f7rvvPm6gbAEGNXINbptNEWIdYqApCN4AAK/64osvdM0112jevHl68MEHtWfPHnXt2tXbZQU0dzbAsVqtMhUXa+uGDcx2A01E8AYAeIXNZtPrr7+uxMREHTt2TOvWrdO8efMUEcGeiJ5glL1Xu6EZ7PCQEM179lml3XSTcnJyPFAZELi4uRIA0CRWSWWVn22yr/HcSlK4nM/uFBQU6O6779bq1as1bNgwLVq0SB07dvRMwaijoX/DwsJC9ezZU1arVf/617903nnnea9YwI8x4w0AaBSL7CtknJNUKqm88lh55e/PyfEKGtu2bVP37t21bt06zZkzR9nZ2YRuL2sle+tJlOyz4FGVv68KCbGxsVq6dKl++OEH3XfffWLODmgcgjcAwG2lcm1ZuqpwXiqpvLxcM2bM0JAhQ9S6dWvt3r1bjzzyiFq14qnIH1xzzTWaNWuW3n33XS1atMjb5QB+iVYTAIBbSiWZ3bzGZrXqL3/+s56cPl133XWX5s2bp9atWQ3a31RUVOi3v/2t9u7dq/379+vXv/61t0sC/ArBGwDgsqoZ7NqWLFqk+XPm6JvDhxXTpo1GpqRo5nPP1egFNplM+iQnR9f/9rceqxfN7+jRo7rqqqv0i1/8Qjk5OTIaWesEcBXv7wEAXFbq4NjLc+boqcce06w//1lHCgu1efduff/dd0oZOlRlZWXV50VGRqo/odvvJSQk6M0339T+/fuVkZHh7XIAv8KMNwDAJVbZb5j8ubNnz+qKTp30yptvKiU9vfp4UVGRul98sZ56/nmNnzixxjUxYtYnEEyaNEkLFizQxo0bdf3113u7HMAvELwBAC4xq+6M95aNG3XLyJE6aTbLYKi5IvR9d9yh8rIyvbFsWY3j7mzeAt9VUlKi3r17Kz8/X5988okuuOACb5cE+DwmHQAALrE6OFaQn6+4+Pg6oVuSOnTsqIL8fJfGgf+JjIzUsmXLdObMGd1xxx2yWvmXBRpC8AYAuMTR26Nx8fEqyM+XxVJ3YcETx48rLj7epXHgn7p166Y5c+Zo48aNmjdvnrfLAXwewRsA4JIQB8d6JyXJaDTq/VWrahwvKirSluxsDRwyxKVx4L8mTZqk0aNH67HHHtOBAwe8XQ7g0+jxBgC4xFGPtyTNe+EFvTJnjhYsXqyBQ4bo2NGjmjppkk6eOKHNubl1lpujxzvw5Ofnq3v37oqJidG+fftYox1wghlvAIBLwp0cf2jaND3x7LN6YupUdWnTRr/t21cJXbrovQ8+cLjGs7Nx4L/i4+P1zjvvKC8vTw899JC3ywF8FjPeAACXubJNfH0MkpgLDVwzZszQc889p+XLlystLc3b5QA+h+ANAHCZs50rXdVa9vCNwFReXq4BAwboyy+/1MGDB3XRRRd5uyTAp9BqAgBwmUGStbhYJSaT29dGiNAd6MLCwrR06VJZrVbddtttDle7AYIZwRsA4LKysjKNHTlSM6dNk7WiwuXrImS/qRKB75e//KUWLFigXbt26ZlnnvF2OYBPIXgDAFxis9l0//33a/v27RqYlKSY0NAGZ7CreroJ3cHltttu0/jx4/X0009rx44dkuwbJ5klmWRvVzJV/p5tdxBM6PEGALjkxRdf1P/+7//qiSee0KxZs6qPWyWVVX62yb5OdyvZVy9hdid4nTt3Tj179tSll1+uZatWKSTc+Xo2BtlfnNGKhEBH8AYANGjt2rVKSUlRWlqali1bplatiNRoWN633yrm/PMVERnp0v8ZWpIQ6PjJCQCo14EDBzRu3Dj17t1bixYtInTDJaWS2v/iF4pq3brG/5kVS5dqUGKiEqKjdVnHjkodPly5O3dKcr5JE/6Llh3/xk9PAIBTx44d06hRoxQXF6e1a9cqMjLS2yXBD1hkD4O1vZKZqcenTNEjM2Yo7+RJfXbkiO6aNEkb1q6tPsespq0VH6iqlvI8J/uLk/LKY+WVvz+npq+zj5ZHqwkAwCGTyaSBAwfq0KFD2rVrl7p37+7tkuAnHAXAwsJCdU1I0PyFCzWmgc112GipplI5fiHjTDC07PjrvSXcxwAAqMNqter222/Xvn37tHbtWkI3XGaV41nXvbm5MpvNGpmS0uAYlspxfDlAeYq7oVs/Oz8Qw7dF9r8TZzP7pfLtm3X5Pw0AqOOJJ57QypUr9eKLL2rUqFHeLgd+pMzJ8dMFBYqLj5fB4FoccjZOMHHWsiPV3ysvBWbLTqlca6epasvxxfsFfPHFAADAi9566y09++yzuueee/Twww97uxz4GWc3+bWLi1NBfr4sFotL4Xvdhg165YUX1Lp1a4cfUVFRLn8tPDxcISEhzftAPcBZcHwlM1MvzZ6tzNde05Drr1d4eLi2bNyoDWvXKql//xrXB0rQC5SZf3q8AQDVduzYoSFDhmjAgAHauHGjwsLCvF0S/EhhYaF+Ki1V2wsucPi1Kzp10oLFi3VTamqDY+3bvVtPTZum4uLi6g+TyaTi4mKVlro3lxkaGupSYHc30FcdCw0NdaseV1hlv2GyNnd65SUpRv7f3lA1g13biqVLNT8zU199+aWiY2LUrUcPPZqRUePFh2S/X8BXXoAQvAEAkqTDhw+rb9++iouL0+7du9W2bVtvlwQfZrPZlJeXp9zcXOXk5Cg3N1eff/65Xlu8WLeMH+/wmpfnzNFfXnhBc19/XYOHDVNYWJi2b9miHdu2adYLL9Q4N0xSlJM/22KxVIfw2qHc0Ye7X7Na3VucLyIiwu3A3tDXouPiFBoVJdWaqd+ycaNuGTlSJ81ml945MMp+s6U/c9Re4mzWP+fDD/X0n/9c41xfulmX4A0A0JkzZ5ScnKyTJ09q9+7duvTSS71dEnxMUVGR9uzZo9zc3OqP06dPS5JiY2OVlJSkpKQkpfzP/+jCSy+tExirLF+yRK/Onau8Q4cUHROjHr166dGMDPVNTq5xnrcCo81mU2lpaYsE+uLiYpnNrjVMvP7WWw5fwCxfskT/79FHlXfihEvjHDl8WIf27VNkZKSioqIUGRlZ4+Pnx1pi5r6pHM38uzvrL/nOzL+vzLwDALzEYrEoPT1d//nPf7R582ZCN2Sz2fT111/XmM3+5JNPqmeCr7jiCo0ZM0bJyclKSkrS5ZdfXr1JjrMWiSrpt92m9Ntua7AG5xvMt6yQkBBFREQoIiJCcXFxzT5+RUWFSkpKnAb2qtDet18/h9e72yv/+eef69ZbbnGptrCwMJcCenMci3RxN1NHN9m6s0LOz8fxhZl/gjcABDGbzaYHH3xQmzdv1ptvvqmBAwd6uyR4gclk0scff1xjNvvHH3+UJEVHR6tv377KyMhQUlKS+vbtq3bt2jkdq5Xs4aIpK2oY5Buzky0hNDRU0dHRio6Orvc8k+yb49TWOylJRqNR69escalX/rrrrtNnn30mk8mkkpKSGh+Ojjk7XlhYqBMnTjg8t7GMRmODAf3+Rx9V31o92+6ukCP5zs6eBG8ACGIvv/yyFixYoGnTpunOO+/0djnwAJvNpiNHjtSYzf7Xv/4li8UelS+99FLdcMMNSkpKUnJysq688kq3WxCMalrw9qVVKLzF2QuP2NhYPT5rlqZOnqxQg6HBXvnYmBi1v/LKFquzqjWnsYHe2bk//fSTjh07ptKyunPe7s76S/ZNdnwBPd4AEKQ2bNigUaNGafTo0Vq5cqVLb/vC/5SWlmr//v3VITsnJ0fHjx+XJEVFRal3797VLSPXXHONzj///Ob5c+X+8m9ScOy66IqGWnZc7ZX3ld7mxnI08+/uCjlS/TfrehLBGwCC0Keffqp+/frpkksu0Y4dO9S6ta/c84+mOnr0aHW7SE5Ojvbv36+yylnDiy++uPomyOTkZF111VVuvV3vLrY6bxpXNoupjy+t5tFYZjlez9ydFXIk31ndheANAEHm5MmT6tu3r8rLy7Vnzx4lJCR4uyQ0UllZmQ4ePFg9m52bm6sjR45IsvfPJiYmVs9mJyUlqUOHDh6vsaEtviXf3uLbm5ytX+0qX1q/urHqm/l3ddZf8p2Zf4I3AAQRs9ms6667TgcPHtSOHTvUq1cvb5fk86yyr4hglb1PNET2J/Bwef6J/OTJkzVmsz/++OPq5em6dOlSYza7R48eCg/31togdfnS36M/oWUnsGb+/f2FEADARTabTRMnTtTu3bu1cuVKQncDGpqprdqOu6Vmai0Wiz799NMas9lff/21JPuyb1dffbXuv//+6rDduXPnFqii+bSSb7zV72+qwnMwt+wE0s26BG8ACBKzZs3SsmXL9Nxzz2ns2LHeLsenuTrLaKn8aI6gk5+fr927d1fPZu/du1fFxfZGgw4dOig5OVn333+/kpOTdfXVVysighgbLIySQhW8LTsG2b/HGjvz70t/H7SaAEAQWLZsmcaNG6c77rhDCxcuVIiTXQXhmbf2Kyoq9MUXX9SYzc7Ly5NkX+e5R48eNXqzL7roIv7NICm4W3YC4WZdgjcABLjdu3dr0KBB6tOnjzZv3iyj0deeinyHo5vZlixapPlz5uibw4cV06aNRqakaOZzz+m8886rc72zm9nOnDlTYzb7o48+0rlz9lvGzj///Bq92YmJiYqK8oWFzwDf4+836xK8ASCAfffdd+rTp4+io6P10UcfKT4+3tsl+bTaN3FVLVm2YPFiDRwyRMeOHtXUSZOUf+qUNu3aVefmRYOkSKtV//73v2vMZn/xxReSpFatWqlbt241ZrN/9atfMZsNuMlfZ/4J3gAQoM6ePat+/frp+++/V25urq644gpvl+TTai9bdvbsWV3RqZNeefNNpaSnVx8vKipS94sv1lPPP6/xEyfWGKOsrEzXdO2qrw8fliS1bdu2xmx27969FRMT44FHA8AXEbwBIABVVFRo9OjR2rRpk7KzszV06FBvl+Tzam/UsWXjRt0ycqROms11Npm57447VF5WpjeWLas5RkmJtqxbp7KiIiUlJenXv/41O4ICqOaL7S8AgCaaOnWqNmzYoAULFhC6XWSt9fuC/HzFxcc73NmxQ8eO+te+fXWOR0RGKiUtzSe2pgbgewjeABBgXnvtNb300kuaMmWK7rvvPm+X41OKi4v1ww8/OPyYNHWqkgcOrD43Lj5eBfn5slgsdcL3iePHFeekX563kQE4Q/AGAD/hys1Emzdv1h/+8AeNGDFCL774oncK9QKbzaazZ886DdVVH2fOnKlzbVxcnDp37ixrRUWN472TkmQ0GvX+qlV1ery3ZGdr5rPPOqyF2yQBOEOPNwD4OFeXzzr+zTdK7NlTF154oXbt2hUwN/HZbDadPn26wVBdVFRU59r27durc+fOTj8SEhIUGRkpqW6PtyTNe+EFvTJnTp1VTU6eOKHNubkOl2Y0ih0aAThG8AYAH+bqhhE2m03mkhI99+STeuQPf9BFF13U0qU1C6vVqvz8fIdB+vvvv6/+tdlc82+hVatW6tixY72hulOnTnWW+6u3FtVc1aTKW3/f+9kAACAASURBVG+8oQVz51av4z1izBg9NXu2zmvb1uE4MfLt5cwAeA/BGwB8VGN2ULRWVCgqNNQndmurqKjQyZMn652lPnr0qMrKympcZzAYlJCQUG+o7tChg8ObHpuq9jre7jLIvokOADhC8AYAH/TzHRRXLF2q+ZmZ+urLLxUdE6NuPXro0YwMJfXv7/R6ZzsoNpfy8nIdP3683lB97NgxVdTqmzYajfUG6s6dO+uCCy7w2hJ8jnaudEdL/70D8G/8fAAAH1TVa/xKZqZemj1bma+9piHXX6/w8HBt2bhRG9aurTd4l6rxP+BLS0t19OjRekP1iRMnVHveJioqSl26dFHnzp01ePBgh6E6Li7Op3dpNMjen+3uOw2qvI4nVQD1YcYbAHxMVa9xYWGhuiYkaP7ChRqTllbnvPsnTFBC5876f888I0nasX277v3d7/TFDz9IctxrbDKZGrxJ8dSpU3X+rDZt2jidoa4K27GxsT4dqt3hbptPhOQT7T0AfBsvzgHAx1R1PO/NzZXZbNbIlBS3x7BYLNqwbZveW7GiRqj+6aef6pzbrl276hDdu3dvhyt/tGnTpomPyr8YJYXKtdVkjOLJFIBr+FkBAD6magfF0wUFTndObIjBYNDR48e1du1ade7cWRdffLEGDBjgMFRHRbHPoiOGyg9X1k8HAFcQvAHAi2w2m3766acaS+glDhyoSy6/XO3i4pzunOiKcbfdpntvv70Fqg4urcS63ACaB8EbAFpI1cYvP1+Puvb61D/88INMJlON6/769tu65PLLq3dOXL9mjW5KTa0zfuvWrWtc++OJEzW+Hhoa2jIPDADQKNxcCaBRgv3td5vNVmfjl58H6qpf1974JTQ0VJ06dapzY+LPf92uY0eVV4bml+fM0V9eeEFzX39dg4cNU1hYmLZv2aId27bpV5deqlfmzNE/cnJUVlam28aM0bEffqi+uZIdFAHAtxC8AbjF1e3L/fmGM5vNplOnTtUbqH/44QeVltbcYDw0NLR64xdHgbpq45eGZqJr76C4fMkSvTp3rvIOHVJ0TIx69OqlRzMy1P3qq3X/HXdoS3a2uvziF7rtzjs1f86celc1AQB4D8EbgMsCYYk1q9WqU6dONdj+4Ww3RWeBunPnzmrfvn2ztXewgyIABB6CNwCXOArdSxYt0vw5c/TN4cOKadNGI1NSNPO553TeeedVn+PJ8G21WmtsUe4oXB89elTl5eU1rgsLC3O4LvXPf+3p3RTZQREAAg/BG0CDHIXAqt7jBYsXa+CQITp29KimTpqk/FOntGnXLoWHh1ef2xwhsKKiojpUO5utPnr0qCyWmvPE4eHh9Qbqzp076/zzz/faFuX1OV1UpLJWrRTp5nJ/vvhOAwCA4A3ABbXbHs6ePasrOnXSK2++qZT09OrjRUVF6n7xxXrq+ec1fuLE6uMNtT1UVFToxIkT9bZ/HDt2rE6oNhqN9QbqqlDtj7splpWVaejQobqsWzc9P2+eWrnYwkLoBgDfxTuRAOplVd1e4z05OTKbzRo1dmyN49HR0Rp6443avnlzjeBdbrNpz8cf67tvv3UYro8fP66KiooaY0VERFSH6IEDB9YJ1F26dFFcXJxfhuqG2Gw23Xvvvfrwww917733KiY0NOBvaAWAYMDPaAD1KnNwrCA/3+mOih06dtS/9u2rcaykpETvvPuuXp4zR5IUGRlZHaSHDBnicLa6Xbt2ARmqXfH8889r0aJFevLJJzVu3DhJ7KAIAIGA4A2gXlYHx+Li453uqHji+HHFxcfXOBYVFaUHpkzR72+/XZ07d1bbtm2DNlQ3ZOXKlXr88cd166236sknn6zxNXZQBAD/xiQJAKesVqsKz56tc7xqR8X3V62qcbyoqEhbsrM1cMiQOtd06txZV111VVDPZDfk448/1vjx45WUlKQ333yTvycACDAEbwDVSkpK9M9//lPPPvusbrzxRsXFxen9tWvrnBcbG6vHnnxS0x54QFs2blR5ebm++/Zb3Zmerk6dO+uW8ePrXEOErN8PP/yg0aNH64ILLtDq1asVEcHcNgAEGlpNgCB26tQp7dq1Szt37tTOnTu1f//+6jWuu3btqrS0NF3ZtatsVqtCai2399C0aWobF6cnpk6tXsd7xJgx+tuSJTIa666rwat854qKijRq1CgVFRUpJydH7du393ZJAIAWwHKCQJCw2WzKy8urDtq7du1SXl6eJPta13369FG/fv3Uv39/JSUlKS4uTlLd7csbi+3LHauoqFBKSorWr1+v9evX64YbbvB2SQCAFuIXM97cyQ+4r6ysTPv3768O2bt27dKpU6ckSe3atVO/fv101113qX///urVq5fDWWrJ/j1mUNO3L+d71bFp06bp/fff1yuvvELoBoAA59Mz3haJtWsBF/3000/Kzc2tDtp79uyR2Wzf5P2SSy6pns3u16+fLrvsMrd2amT78pbx17/+Vffee68eeOAB/eUvf/F2OQCAFuazwbtUktmN89mtDcHEZrPpu+++q+7N3rVrlz777DNJksFgUM+ePatDdr9+/dShQ4cm/5nufk9W4XvTsQ8++EA33HCDhg4dqvfee8/hmugAgMDik8GbJ/j/os0GkmSxWPTJJ59Uh+ydO3fq2LFjkqSYmBglJydXB+0+ffqodev6NmhvPF4QN48vv/xS11xzjbp06aJdu3apTZs23i4JAOABPhe8nb2lvWLpUs3PzNRXX36p6JgYdevRQ49mZCipf/8a5wXKW9q02QS3oqIi7d69uzpo7969W0VFRZKkLl26VIfs/v376ze/+Y1CQ0M9Vhv/N5smPz9fffv2VVFRkfbs2aOLLrrI2yUBADzE554XSx0ceyUzUy/Nnq3M117TkOuvV3h4uLZs3KgNa9fWCd6l8sEH5SZXZxUtlR/MKvq/Y8eO1ZjNPnjwoCoqKhQSEqKrrrpKt99+e3XYvvDCC71aq0FsX95YpaWlGjt2rI4ePart27cTugEgyPjUjLejZcsKCwvVNSFB8xcu1Ji0NJfG8edly2izaTm+EhStVqu++OKLGkH722+/lSRFRkbqmmuuqZ7NvuaaaxQbG+vB6tBSbDabJkyYoLfeekvLli3T//zP/3i7JACAh/nU5HCZg2N7c3NlNps1MiXFrXH8cc83i5yH7oZabcySQuVj/6A+oqHWiKp3SVqqNaKkpER79+6tDtk5OTk6c+aMJKl9+/bq16+fHnzwQfXv3189evRQWFhYC1QBb3vuuef01ltv6Y9//COhGwCClE/lNKuDY6cLChQXH+/WHf+OxvEHjtpsJNdbbQKhzaa5eaNtp2o3yKqgvW/fvurdIC+//HKlpqZWt4386le/UkgIm6kHuhUrVigjI0O33XabnnjiCW+XAwDwEp9qNSlW3VnJLRs36paRI3XSbHY5fG9av173jR+vyMhIhx9RUVFOv9aYc91ZD9kZZ7sDuttq489tNs3NE207NptNX331VY3dIP/9739Lsu8GmZiYWB2yk5OTFR8f34iK4M/27NmjgQMH6uqrr9YHH3ygiAh/fD8OANAcfCp4mySV1zpWWFioKzp10oLFi3VTaqpL43y6b5+WL14sk8mkkpKSBj9MJlP1RiONER4e3uRA/5vERF3RvbtCa724cPeFh1H+2WbT3Orb8MWVFXKcrY5TVlamAwcO1Fg/u2o3yLZt29bYpCYxMZGQFeSOHDmiPn36KCoqSh999JHOP/98b5cEAPAin+pMcDRTGxsbq8dnzdLUyZMVajBo8LBhCgsL0/YtW7Rj2zbNeuGFOtck9uql/r16ufVn22w2mc3mesO5KwG+9rHCwkKdOHHC4Xk/9/pbb+k3Dmp2t9Vm67ZtWvL3v6t169aKiopSVFSUw1/X9/WIiAi/b39orradM2fO1NgN8qOPPqp+kfbLX/5Sw4cPrw7bl19+ebO8+4HAcO7cOY0aNUolJSXaunUroRsA4Fsz3s7aLSRp+ZIlenXuXOUdOqTomBj16NVLj2ZkqG9ycp1z/aHdwmazqbS0tDqMG2JjZXSw6Ym7M967tm/Xg7//vYqLi2UymVRcXKyKigq362sonDf11+Hh4S0W7pujbae8vFxjf/tb7dyxQzabTaGhoerZs2eNGe2OHTu2SP3wfxUVFbrpppu0ceNGbdiwQcOGDfN2SQAAH+BTwVty3OftDoPsbQL+xlGbjeR+q02YpKhax8rKymQymaqDeEv82t3/RqGhodVh3NWZeJd/fd55UkSEVCvYu/MipsRkUtaSJTp9/Lj69++vPn36KDo62q3HiOD18MMP66WXXtKrr76q+++/39vlAAB8hE+1mkj2HuWmBG9/Xcva2Qy9u602jsYJDw9XeHi4zjvvvGavW/rv7L2rQb2h806cOOHwGle9/tZbumX8+DrH3WnbiYyK0sS7767zIgZoyGuvvaaXXnpJDz30EKEbAFCDzwVvg+w3BzZ2NQqfe0AuCpfzvuQHHn1U7Tt00IvPPKN7brutRquNo3E8LSQkRBEREYqIiFC7du1a5M+wWq0ym80uhfskB+1HktQuLk4F+fmyWCwuhW+feisIfmHz5s36wx/+oBtvvFFz5szxdjkAAB/jc60mVdxdCi4Qdm4M1jab5taSbTuAM4cOHVJSUpIuvPBC7dq1SzExMd4uCQDgY3z2HkSjnC/p9nNVYdPfQ7fU9McQCH8HzcGVtp11a9bIZDKpvLxcm7OzNXPaNJfHAWo7deqURowYoYiICK1bt47QDQBwyGdnvH/OKvs28FbZ3/4PkT0UhSvwwpEnNn0JdPWtjiO5vkKOP6yOA+8rLS3VkCFDtG/fPv3zn/9Unz59vF0SAMBH+UXwDjbB2GbT3GjbgSfYbDbdfvvteuedd/Tuu+8qPT3d2yUBAHyYv96LGNCMkkJlD+D1hUdD5bn8I9YVrKvjwLP+9Kc/6Z133tHTTz9N6AYANIgZbx8XTG02za0xbTsWi0XRBgPBGw1+7y1fvly33HKLxo8fr8WLF/v9bq8AgJZH8EZAczl822wyl5bqyWnTNPnuu9WtW7cWrgy+yqKG3206d/q00kaPliEkRFu2bJHRyEs1AEDDCN4IeK4EKYOkovx89bzqKkVHR2vv3r2KjY31TIHwGa6+UKuoqFBZaanCKioUxwomAAAX0a2AgFd1o2SM7L3bYZXHwip/H1P59fbx8Xr33Xf19ddfa+LEieI1aXBxpzUpNDRUkVFRMsTEON34CgCA2gjeCBqtZF8BJkr2oB1V+fuffxMMGDBAs2fP1qpVqzR37lwvVAlvsKhu6F6yaJGSu3VTx6go/bpDBz1y//06c+ZMnWvNatqNvACA4EGrCVCLzWbTzTffrPfee0/bt29X//79vV0SWljt5SdfnjNHf3nhBS1YvFgDhwzRsaNHNXXSJOWfOqVNu3YpPDy8xvUsPwkAcAXBG3CgsLBQiYmJMplMOnDggC644AJvl4QWUnvDpbNnz+qKTp30yptvKuVnSwQWFRWp+8UX66nnn9f4iRPrjMOGSwCAhvA8ATgQGxurrKwsnT59WrfeeqsqKiq8XRJaSFmt3+/JyZHZbNaosWNrHI+OjtbQG2/U9s2bXRoHAIDaCN6AE927d9eCBQu0detWzZw509vloIVYa/2+ID9fcfHxMhjqbk3VoWNHFeTnuzQOAAC1EbyBekyYMEG///3v9eyzz2rdunXeLgctoHavXVx8vAry82Wx1L1l8sTx44qLj3dpHAAAaiN4Aw14+eWX1bNnT40fP17ffPONt8tBM6u932TvpCQZjUa9v2pVjeNFRUXakp2tgUOGuDQOAAC1EbyBBkRERCgrK0s2m02pqakym93diB6+rPYPwdjYWD325JOa9sAD2rJxo8rLy/Xdt9/qzvR0dercWbeMH+/SOAAA1MZzBeCCX/7yl3rrrbe0f/9+PfTQQ94uB80oXKqzWdJD06bpiWef1RNTp6pLmzb6bd++SujSRe998IHT7eHDHR4FAOC/WE4QcMP06dP1/PPPa/Hixbr99tu9XQ6awY4dO/RjUZGuGzZMoaGhjRqDdbwBAK5gxhtwwzPPPKNBgwbpvvvu06effurtctAEFotFM2fO1KBBg/TW3/7WpB5tx3PgAADUxIw34KYTJ07o6quvVnR0tD7++GO1adPG2yXBTd98841uu+025ebm6o477tDLL7+s8JiYOtvGuyJCBG8AgGuY8Qbc1KFDB7377rv6+uuvNXHixDr9wfBty5YtU48ePfT5559r6dKlWrRokWJiYmSUPUS7g9ANAHAHwRtohAEDBmj27NlauXKlXnrpJW+XAxecO3dOEyZM0Lhx43TllVfqX//6l2699dYa5xhl79Wuu3VOTVU93YRuAIA7aDUBGslms+nmm2/W+++/r+3bt6tfv37eLglO7N27V+PGjdPXX3+tjIwMzZw50+HOlD9nlX0beKvsm+OEyD5TES5mLAAAjUPwBpqgsLBQiYmJMplMOnDggC644AICmw+xWq168cUXlZGRoQ4dOmjJkiW69tprvV0WACBIkQOAJoiNjVVWVpZOnz6tJ556SkVWq85JKpVULslS+blU0jlJxZXH0PKOHTumoUOH6rHHHtNNN92kTz75hNANAPAqZryBZrA9J0eXde+uiMhItWrV8OtZbsprWe+9954mTpyokpISzZs3T3fddZdCQtjUHQDgXQRvoIlKJZah8xElJSWaOnWqXn31VfXo0UPLli3T5Zdf7u2yAACQRKsJ0CQWOQ/dK5Yu1aDERCVER+uyjh2VOny4cnfurP66WbSdNKfPPvtMffr00auvvqqHH35Yu3fvJnQDAHwKwRtoglInx1/JzNTjU6bokRkzlHfypD47ckR3TZqkDWvXunQ9XGez2TR//nwlJibqxx9/VHZ2tjIzM2U08n4CAMC30GoCNJJV9hsmayssLFTXhATNX7hQY9LSGhwnRrwCbqz8/HxNnDhR77//voYPH66FCxeqffv23i4LAACHeL4HGqnMyfG9ubkym80amZLSpHFQvy1btuiqq67Spk2bNHfuXK1bt47QDQDwaQRvoJGsTo6fLihQXHx8gxu0NDQOHCsrK9Njjz2mYcOGKTY2Vh999JGmTJni0moyAAB4k2vJAEAdznq02sXFqSA/XxaLxaXwTa+X67766iuNGzdOH3/8se655x5lZmaqdevW3i4LAACXMEUENJKzVaF7JyXJaDRq/Zo1Lo1TdM5Rpzh+zmazafHixerZs6cOHz6srKwsvf7664RuAIBfIXgDjeTsmyc2NlaPz5qlqZMna92aNTKZTCovL9fm7GzNnDatxrkmk0l/+uMf1atXLz333HP6z3/+0/KF+5kzZ85o3LhxmjBhghITE3Xw4EHdfPPN3i4LAAC3saoJ0EjOVjWpsnzJEr06d67yDh1SdEyMevTqpUczMtQ3Obn6HJvVqncWLNA7b7+tjz76SJLUvXt3paamKi0tTZdddlnLPggfl5OTo3HjxumHH37QH//4R02fPl2hoaHeLgsAgEYheANNUKymbYJjkFTVLHHkyBGtXLlSWVlZysnJkST95je/qQ7hXbt2bWK1/qOiokJ/+tOfNGvWLHXp0kVLly5VUlKSt8sCAKBJCN5AE1hkD9+N1VqO73A+evRodQjfuXOnbDabrrjiCqWlpSk1NVW/+c1vFBLirMvcvx05ckS/+93vtGPHDo0bN06vvvqqYmNjvV0WAABNRvAGmqhUzreNr0+EJFf2Vjx+/LhWrVqlrKwsffjhh7JarbrsssuUmpqq1NRUde/ePWBCeFZWlu6++25ZLBbNnz9f48ePD5jHBgAAwRtoBu6Gb1dDd20nT57U6tWrlZWVpW3btslqteqSSy6pDuFXX321XwbV4uJiTZkyRX//+9/Vu3dvLV26VJdccom3ywIAoFkRvIFmYpE9gNfX822QPXA3xwL6p06d0po1a5SVlaUPPvhAFRUVuvjii6tDeO/evf0ihB84cEC33nqr8vLy9Nhjj2nWrFkKCwvzdlkAADQ7gjfQzKyybwNvlX1znBDZlx4MV8ut31lQUKC1a9dqxYoV2rJliywWiy688MLqGzP79Onjczs7Wq1WvfTSS5o+fbrOP/98vf322xo8eLC3ywIAoMUQvIEA89NPP2nt2rXKysrSP/7xD5WXl6tz5866+eablZaWpqSkpGYN4Y15oXHixAlNmDBBmzZt0ujRo/XGG28oPj6+2WoCAMAXEbyBAHbmzBm9//77ysrK0qZNm1RaWqqOHTtWh/B+/fo1el3sxrbWZGdna8KECTp79qwyMzN13333+UVLDAAATUXwBoLE2bNntW7dOmVlZSk7O1tms1nt27fX2LFjlZaWpgEDBshgcK37vDE3k6q0VI899pjmzZunbt26admyZbryyisb8UgAAPBPBG8gCJ07d04bNmxQVlaW1q9fr5KSEp1//vkaO3asUlNTNWjQIKchvDHLJ1qtVs199lk9/cQTeuCBB/TCCy8oIiKiyY8DAAB/QvAGglxxcbGys7OVlZWldevWqbi4WHFxcUpJSVFqaqoGDx5cvcpIfRsGrVi6VPMzM/XVl18qOiZG3Xr00KMZGUrq31+SZDKZ9O99+zRowADPPDAAAHwMwRtANZPJpE2bNmnFihV6//33VVRUpLZt22rMmDFKTU3VgOuvl9VBT/grmZl6afZsZb72moZcf73Cw8O1ZeNG5Xz4oZ7+858lSTabTWEhIWrt6QcFAICPIHgDcMhsNusf//iHVqxYoffee0/hRqM+O3KkTotIYWGhuiYkaP7ChRqTltbguDFquWUVAQDwZQRvAA0qLS3Voa+/VsKvfqXw8PAaX9uycaNuGTlSJ81ml27ONKryZksAAIJMc2ygByDAGY1G/fqKK1Tu4GunCwoUFx/v8ooo1uYtDQAAv8E7vgBc4uytsXZxcSrIz5fFUt+K3g2PAwBAoCN4A3CJsy1ueiclyWg0av2aNU0aBwCAQEfwBuASZz8sYmNj9fisWZo6ebLWrVkjk8mk8vJybc7O1sxp01weBwCAQMfNlQBcYpV0rp6vL1+yRK/Onau8Q4cUHROjHr166dGMDPVNTq5xHquaAACCFcEbgMuKZd9Ep7EMEut4AwCCFhNPAFxm9PL1AAD4M4I3AJcZ1Pg1uCPE+qUAgOBG8AbglsZsgBMhZrsBAKDHG0CjWCSVSjKVlclaUaGIyMg65xhkD9zMdAMAwIw3gEYySDJaLOpz+eXasm6dwiqPhcketmNkv5GS0A0AgB3PiQAa7cCBA/r2m28UWlGhKG8XAwCAj2PGG0Cjbd26VZJ03XXXebkSAAB8H8EbQKNt3bpVV155pdq3b+/tUgAA8HkEbwCNUlZWph07dmjw4MHeLgUAAL9A8AbQKB999JFKSkoI3gAAuIjgDaBRtm7dqpCQEA0cONDbpQAA4BdYxxtAowwaNEhFRUX6+OOPvV0KAAB+gRlvAG4zmUzKzc2lzQQAADcQvAG4LScnR2VlZQRvAADcQPAG4LatW7fKYDCof//+3i4FAAC/QY83ALddc801Cg0N1a5du7xdCgAAfoMZbwBuKSws1N69e2kzAQDATQRvAG7ZsWOHrFYrwRsAADcRvAG4ZevWrTIajUpKSvJ2KQAA+BWCNwC3bN26Vf369VNERIS3SwEAwK8QvAG4LD8/XwcPHqTNBACARiB4A3DZ9u3bJYngDQBAIxC8Abhs69atio6OVmJiordLAQDA7xC8Abhs69atuvbaaxUWFubtUgAA8DsEbwAuOXbsmP7973/TZgIAQCMRvAG4ZNu2bZLo7wYAoLEI3gBcsnXrVrVt21bdu3f3dikAAPglgjcAl2zdulXXXXedWrXixwYAAI3BMyiABn3zzTf69ttvdd1113m7FAAA/BbBG0CDtm7dKon+bgAAmsLg7QIA+B6rpLLKzzZJXS67TDP++EdddsUV3i0MAAA/FmKz2WzeLgKAb7BIKq38XFtZaanCjUYZJBnFq3YAANxF8AYgyR64zW6cHyF7AAcAAK6hxxuA26FbleeXtkAtAAAEKoI3EOQschy6VyxdqkGJiUqIjtZlHTsqdfhw5e7cWeMcsxy3pQAAgLoI3kCQczRr/Upmph6fMkWPzJihvJMn9dmRI7pr0iRtWLvWpesBAEBd9HgDQcwq6VytY4WFheqakKD5CxdqTFqaS+PEiFfxAAA0hOdKIIiVOTi2NzdXZrNZI1NSmjQOAACoieANBDGrg2OnCwoUFx8vg8H1BQMdjQMAAGoieANBzFGfWbu4OBXk58ticf22SfrVAABoGMEbCGIhDo71TkqS0WjU+jVrmjQOAACoieANBDFHPwBiY2P1+KxZmjp5statWSOTyaTy8nJtzs7WzGnTXB4HAADUxKomQBBztKpJleVLlujVuXOVd+iQomNi1KNXLz2akaG+ycl1zmVVEwAAGkbwBoJcsZq2CY5BUutmqgUAgEDGJBUQ5Ixevh4AgGBB8AaCnEFSRCOvjai8HgAANIzgDUBG/Sx8u9h9FiFmuwEAcAfBG4Ake4huLenU8eMyl5SooqLC4XlVPd2EbgAA3EPwBlDNIGnO00+rT9euCrfZFFZ5LEz2oB0je+imvQQAAPfx/Amgms1m04YNG3T11Vcr2o0t4wEAQMOY8QZQ7dChQzpy5IhuvPFGb5cCAEDAIXgDqLZhwwZJ0vDhw71cCQAAgYcNdABUGzx4sPLz8/XJJ594uxQAAAIOM94AJElnz57Vjh07aDMBAKCFELwBSJK2bNkii8VC8AYAoIUQvAFIkrKzs9WmTRslJSV5uxQAAAISwRtA9TKCw4YNU1hYmLfLAQAgIBG8AeiTTz7RsWPHaDMBAKAFEbwBVC8jeMMNN3i5EgAAAhfLCQLQtddeq+LiYu3bt8/bpQAAELCY8QaC3E8/915lwQAAArNJREFU/aScnBw2zQEAoIURvIEgt3nzZlVUVNDfDQBACyN4A0Fuw4YNateunfr27evtUgAACGgGbxcAwHOsksoqP9skyWZTwsUX6+bUVIWGhnq1NgAAAh03VwJBwCKptPJzbSaTScbwcBkNBhnFq3EAAFoKwRsIcKWSzG6cHyHJ2EK1AAAQzOjxBgKYu6FbleeXtkAtAAAEO4I3EKAsch66VyxdqkGJiUqIjtZlHTsqdfhw5e7cWf11sxy3pQAAgMYjeAMBytms9SuZmXp8yhQ9MmOG8k6e1GdHjuiuSZO0Ye1al64HAACNQ483EICsks45OF5YWKiuCQmav3ChxqSlNThOjHh1DgBAc+E5FQhAZU6O783Nldls1siUlCaNAwAA3EfwBgKQ1cnx0wUFiouPl8Hg2qKBzsYBAADuI3gDAchZ/1i7uDgV5OfLYnHt1kn60AAAaD4EbyAAhTg53jspSUajUevXrGnSOAAAwH0EbyAAOfvGjo2N1eOzZmnq5Mlat2aNTCaTysvLtTk7WzOnTXN5HAAA4D5WNQECkLNVTaosX7JEr86dq7xDhxQdE6MevXrp0YwM9U1OrnEeq5oAANB8CN5AgCpW0zbBMUhq3Uy1AAAAJrOAgGX08vUAAKAmgjcQoAySIhp5bUTl9QAAoPkQvIEAZpT74TtCzHYDANAS6PEGgoBFUqnq7/k2yB64mekGAKBlELyBIGKVfRt4q+yb44TI/rZXuHj7CwCAlkbwBgAAADyASS4AAADAAwjeAAAAgAcQvAEAAAAPIHgDAAAAHkDwBgAAADyA4A0AAAB4AMEbAAAA8ACCNwAAAOABBG8AAADAAwjeAAAAgAcQvAEAAAAPIHgDAAAAHkDwBgAAADyA4A0AAAB4AMEbAAAA8ACCNwAAAOABBG8AAADAA/4/LqhWNsXQE00AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "i = 0\n",
    "mol = to_molecule(dataset[i])\n",
    "plt.figure(figsize=(10, 5))\n",
    "draw_molecule(mol)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50c1eb31",
   "metadata": {},
   "source": [
    "## Train and evaluate a GNN classification model\n",
    "\n",
    "As our classifier we use a variation of a [Graph Isomorphism Network](https://arxiv.org/abs/1810.00826) incorporating edge (bond) as well as node (atom) features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "4aeb1f79",
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch.nn as nn\n",
    "import torch.nn.functional as F\n",
    "from torch_geometric.data.batch import Batch\n",
    "from torch_geometric.nn import MessagePassing, global_add_pool, global_max_pool, global_mean_pool, LayerNorm\n",
    "from ogb.graphproppred.mol_encoder import AtomEncoder, BondEncoder\n",
    "\n",
    "\n",
    "class GINConv(MessagePassing):\n",
    "    def __init__(self, emb_dim: int) -> None:\n",
    "        super().__init__(aggr='add')\n",
    "        self.mlp = nn.Sequential(\n",
    "            nn.Linear(emb_dim, 2 * emb_dim),\n",
    "            nn.BatchNorm1d(2 * emb_dim),\n",
    "            nn.ReLU(),\n",
    "            nn.Linear(2 * emb_dim, emb_dim)\n",
    "        )\n",
    "        self.eps = nn.Parameter(torch.Tensor([0.]))\n",
    "        self.bond_encoder = BondEncoder(emb_dim=emb_dim)  # encode edge features\n",
    "\n",
    "    def forward(self, x: torch.Tensor, edge_index: torch.Tensor, edge_attr: torch.Tensor) -> torch.Tensor:\n",
    "        edge_emb = self.bond_encoder(edge_attr)\n",
    "        return self.mlp((1 + self.eps) * x + self.propagate(edge_index, x=x, edge_attr=edge_emb))\n",
    "\n",
    "    def message(self, x_j: torch.Tensor, edge_attr: torch.Tensor) -> torch.Tensor:\n",
    "        return x_j + edge_attr\n",
    "\n",
    "    def update(self, aggr_out: torch.Tensor) -> torch.Tensor:\n",
    "        return aggr_out\n",
    "\n",
    "\n",
    "class GIN(nn.Module):\n",
    "    def __init__(self, n_layer: int = 5, emb_dim: int = 64, n_out: int = 2, dropout: float = .5,\n",
    "                 jk: bool = True, residual: bool = True, pool: str = 'add', norm: str = 'batch') -> None:\n",
    "        super().__init__()\n",
    "        self.n_layer = n_layer\n",
    "        self.jk = jk  # jumping-knowledge\n",
    "        self.residual = residual  # residual/skip connections\n",
    "        self.atom_encoder = AtomEncoder(emb_dim=emb_dim)  # encode node features\n",
    "        self.convs = nn.ModuleList([GINConv(emb_dim) for _ in range(n_layer)])\n",
    "        norm = nn.BatchNorm1d if norm == 'batch' else LayerNorm\n",
    "        self.bns = nn.ModuleList([norm(emb_dim) for _ in range(n_layer)])\n",
    "        if pool == 'mean':\n",
    "            self.pool = global_mean_pool\n",
    "        elif pool == 'add':\n",
    "            self.pool = global_add_pool\n",
    "        elif pool == 'max':\n",
    "            self.pool = global_max_pool\n",
    "        pool_dim = (n_layer + 1) * emb_dim if jk else emb_dim\n",
    "        self.linear = nn.Linear(pool_dim, n_out)\n",
    "        self.dropout = nn.Dropout(p=dropout)\n",
    "\n",
    "    def forward(self, data: Batch) -> torch.Tensor:\n",
    "        x, edge_index, edge_attr, batch = data.x, data.edge_index, data.edge_attr, data.batch\n",
    "        # node embeddings\n",
    "        hs = [self.atom_encoder(x)]\n",
    "        for layer in range(self.n_layer):\n",
    "            h = self.convs[layer](hs[layer], edge_index, edge_attr)\n",
    "            h = self.bns[layer](h)\n",
    "            if layer < self.n_layer - 1:\n",
    "                h = F.relu(h)\n",
    "            if self.residual:\n",
    "                h += hs[layer]\n",
    "            hs += [h]\n",
    "        # graph embedding and prediction\n",
    "        if self.jk:\n",
    "            h = torch.cat([h for h in hs], -1)\n",
    "        h_pool = self.pool(h, batch)\n",
    "        h_drop = self.dropout(h_pool)\n",
    "        h_out = self.linear(h_drop)\n",
    "        return h_out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "64275293",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "device: cuda\n"
     ]
    }
   ],
   "source": [
    "device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\n",
    "print(f'device: {device}')\n",
    "\n",
    "n_layer = 5\n",
    "emb_dim = 300\n",
    "n_out = 1\n",
    "dropout = .5\n",
    "jk = True\n",
    "residual = False\n",
    "pool = 'mean'\n",
    "norm = 'batch'\n",
    "\n",
    "model = GIN(n_layer, emb_dim, n_out, dropout, jk, residual, pool, norm).to(device)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "332f2f6c",
   "metadata": {},
   "source": [
    "Train and evaluate the model. Evaluation is done using [ROC-AUC](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.roc_auc_score.html). If you already have a trained model saved, you can directly load it by specifying the `load_path`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "e9ad612b",
   "metadata": {},
   "outputs": [],
   "source": [
    "load_path = 'gnn'  # set to None if no pretrained model available"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "431050be",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "train ROC-AUC: 0.961 -- loss: 0.070\n",
      "val ROC-AUC: 0.829 -- loss: 0.084\n",
      "test ROC-AUC: 0.753 -- loss: 0.163\n"
     ]
    }
   ],
   "source": [
    "from ogb.graphproppred import Evaluator\n",
    "from tqdm import tqdm\n",
    "\n",
    "criterion = nn.BCEWithLogitsLoss()\n",
    "optim = torch.optim.Adam(model.parameters(), lr=.001)\n",
    "evaluator = Evaluator(name=dataset_name)  # ROC-AUC for ogbg-molhiv\n",
    "\n",
    "\n",
    "def train(loader: DataLoader, verbose: bool = False) -> None:\n",
    "    dl = tqdm(loader, total=len(loader)) if verbose else loader\n",
    "    model.train()\n",
    "    for data in dl:\n",
    "        data = data.to(device)\n",
    "        optim.zero_grad()\n",
    "        y_hat = model(data)\n",
    "        is_labeled = data.y == data.y\n",
    "        loss = criterion(y_hat[is_labeled], data.y[is_labeled].float())\n",
    "        loss.backward()\n",
    "        optim.step()\n",
    "        if verbose:\n",
    "            dl.set_postfix(dict(loss=loss.item()))\n",
    "\n",
    "\n",
    "def evaluate(loader: DataLoader, split: str, verbose: bool = False) -> float:\n",
    "    dl = tqdm(loader, total=len(loader)) if verbose else loader\n",
    "    model.eval()\n",
    "    y_pred, y_true = [], []\n",
    "    for data in dl:\n",
    "        data = data.to(device)\n",
    "        with torch.no_grad():\n",
    "            y_hat = model(data)\n",
    "        y_pred.append(y_hat.cpu())\n",
    "        y_true.append(data.y.float().cpu())\n",
    "    y_true = torch.cat(y_true, dim=0)\n",
    "    y_pred = torch.cat(y_pred, dim=0)\n",
    "    loss = criterion(y_pred, y_true)\n",
    "    input_dict = dict(y_true=y_true, y_pred=y_pred)\n",
    "    result_dict = evaluator.eval(input_dict)\n",
    "    print(f'{split} ROC-AUC: {result_dict[\"rocauc\"]:.3f} -- loss: {loss:.3f}')\n",
    "    return result_dict[\"rocauc\"]\n",
    "\n",
    "\n",
    "if load_path is None or not os.path.isdir(load_path):\n",
    "    epochs = 150\n",
    "    rocauc_best = 0.\n",
    "    save_path = 'gnn'\n",
    "    for epoch in range(epochs):\n",
    "        print(f'\\nEpoch {epoch + 1} / {epochs}')\n",
    "        train(dl_tr)\n",
    "        _ = evaluate(dl_tr, 'train')\n",
    "        rocauc = evaluate(dl_val, 'val')\n",
    "        if rocauc > rocauc_best and os.path.isdir(save_path):\n",
    "            print('Saving new best model.')\n",
    "            rocauc_best = rocauc\n",
    "            torch.save(model.state_dict(), os.path.join(save_path, 'model.dict'))\n",
    "        _ = evaluate(dl_te, 'test')\n",
    "    load_path = save_path\n",
    "\n",
    "\n",
    "# load (best) model\n",
    "model.load_state_dict(torch.load(os.path.join(load_path, 'model.dict')))\n",
    "_ = evaluate(dl_tr, 'train')\n",
    "_ = evaluate(dl_val, 'val')\n",
    "_ = evaluate(dl_te, 'test')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d10f97b7",
   "metadata": {},
   "source": [
    "## Detect drift\n",
    "\n",
    "### Prediction distribution drift\n",
    "\n",
    "We will first detect drift on the prediction distribution of the GIN model. Since the binary classification model returns continuous numerical univariate predictions, we use the [Kolmogorov-Smirnov drift detector](https://docs.seldon.io/projects/alibi-detect/en/stable/cd/methods/ksdrift.html). First we define some utility functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "4b893d04",
   "metadata": {},
   "outputs": [],
   "source": [
    "from torch_geometric.data import Batch, Data\n",
    "from typing import Dict, List, Union\n",
    "\n",
    "\n",
    "labels = ['No!', 'Yes!']\n",
    "def make_predictions(dd, xs: Dict[str, List[Data]]) -> None:\n",
    "    for split, x in xs.items():\n",
    "        preds = dd.predict(x)\n",
    "        dl = DataLoader(x, batch_size=32, shuffle=False)\n",
    "        _ = evaluate(dl, split)\n",
    "        print('Drift? {}'.format(labels[preds['data']['is_drift']]))\n",
    "        if isinstance(preds[\"data\"][\"p_val\"], (list, np.ndarray)):\n",
    "            print(f'p-value: {preds[\"data\"][\"p_val\"]}')\n",
    "        else:\n",
    "            print(f'p-value: {preds[\"data\"][\"p_val\"]:.3f}')\n",
    "        print('')\n",
    "        \n",
    "\n",
    "def sample(split: str, n: int, ) -> List[Data]:\n",
    "    idx = np.random.choice(split_idx[split].numpy(), size=n, replace=False)\n",
    "    return [dataset[i] for i in idx]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff360014",
   "metadata": {},
   "source": [
    "Because we pass lists with `torch_geometric.data.Data` objects to the detector, we need to preprocess the data using the `batch_fn` into `torch_geometric.data.Batch` objects which can be fed to the model. Then we detect drift on the model prediction distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "965f880d",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:alibi_detect.cd.utils:Input shape could not be inferred. If alibi_detect.models.tensorflow.embedding.TransformerEmbedding is used as preprocessing step, a saved detector cannot be reinitialized.\n"
     ]
    }
   ],
   "source": [
    "from alibi_detect.cd import KSDrift\n",
    "from alibi_detect.utils.pytorch import predict_batch\n",
    "from functools import partial\n",
    "\n",
    "\n",
    "def batch_fn(data: Union[List[Data], Batch]) -> Batch:\n",
    "    if isinstance(data, Batch):\n",
    "        return data\n",
    "    else:\n",
    "        return Batch().from_data_list(data)\n",
    "\n",
    "\n",
    "preprocess_fn = partial(predict_batch, model=model, device=device, preprocess_fn=batch_fn, batch_size=32)\n",
    "dd = KSDrift(x_ref, p_val=.05, preprocess_fn=preprocess_fn)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f003b84",
   "metadata": {},
   "source": [
    "Since the dataset is heavily imbalanced, we will test the detectors on a sample which oversamples from the minority class (molecules which inhibit HIV virus replication):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "1a7cf139",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# instances: 630 -- # class 1: 130\n"
     ]
    }
   ],
   "source": [
    "split = 'test'\n",
    "x_imb = sample(split, 500)\n",
    "n = 0\n",
    "for i in split_idx[split]:\n",
    "    if dataset[i].y[0].item() == 1:\n",
    "        x_imb.append(dataset[i])\n",
    "        n += 1\n",
    "print(f'# instances: {len(x_imb)} -- # class 1: {n}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "5e8bfb9f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "H0 ROC-AUC: 0.763 -- loss: 0.208\n",
      "Drift? No!\n",
      "p-value: [0.84936297]\n",
      "\n",
      "test sample ROC-AUC: 0.736 -- loss: 0.150\n",
      "Drift? No!\n",
      "p-value: [0.13626936]\n",
      "\n",
      "imbalanced sample ROC-AUC: 0.749 -- loss: 0.940\n",
      "Drift? Yes!\n",
      "p-value: [4.5092685e-05]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "xs = {'H0': x_h0, 'test sample': sample('test', 500), 'imbalanced sample': x_imb}\n",
    "make_predictions(dd, xs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43a9054e",
   "metadata": {},
   "source": [
    "As expected, prediction distribution shift is detected for the imbalanced sample but not for the random test sample with similar label distribution as the reference data.\n",
    "\n",
    "### Prediction uncertainty drift\n",
    "\n",
    "The [model uncertainty drift detector](https://docs.seldon.io/projects/alibi-detect/en/stable/cd/methods/modeluncdrift.html) can pick up when the model predictions drift into areas of changed uncertainty compared to the reference data. This can be a good proxy for drift which results in model performance degradation. The uncertainty is estimated via a Monte Carlo estimate ([MC dropout](https://arxiv.org/abs/1506.02142)). We use the [RegressorUncertaintyDrift detector](https://docs.seldon.io/projects/alibi-detect/en/stable/cd/methods/modeluncdrift.html) since our binary classification model returns 1D logits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "4e701bff",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:alibi_detect.cd.utils:Input shape could not be inferred. If alibi_detect.models.tensorflow.embedding.TransformerEmbedding is used as preprocessing step, a saved detector cannot be reinitialized.\n"
     ]
    }
   ],
   "source": [
    "from alibi_detect.cd import RegressorUncertaintyDrift\n",
    "\n",
    "dd = RegressorUncertaintyDrift(x_ref, model=model, backend='pytorch', p_val=.05, n_evals=100,\n",
    "                               uncertainty_type='mc_dropout', preprocess_batch_fn=batch_fn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "fbbff246",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "H0 ROC-AUC: 0.763 -- loss: 0.208\n",
      "Drift? No!\n",
      "p-value: [0.52494323]\n",
      "\n",
      "test sample ROC-AUC: 0.736 -- loss: 0.150\n",
      "Drift? Yes!\n",
      "p-value: [0.]\n",
      "\n",
      "imbalanced sample ROC-AUC: 0.749 -- loss: 0.940\n",
      "Drift? Yes!\n",
      "p-value: [0.]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "make_predictions(dd, xs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b765295",
   "metadata": {},
   "source": [
    "Although we didn't pick up drift in the GIN model prediction distribution for the test sample, we can see that the model is less certain about the predictions on the test set, illustrated by the lower ROC-AUC.\n",
    "\n",
    "### Input data drift using the MMD detector\n",
    "\n",
    "We can also more detect drift on the input data by encoding the data with a randomly initialized GNN to extract graph embeddings. Then we apply our detector of choice, e.g. the [MMD detector](https://docs.seldon.io/projects/alibi-detect/en/stable/cd/methods/mmddrift.html) on the extracted embeddings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "897077b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "class Encoder(nn.Module):\n",
    "    def __init__(self, n_layer: int = 1, emb_dim: int = 64, jk: bool = True, \n",
    "                 residual: bool = True, pool: str = 'add', norm: str = 'batch') -> None:\n",
    "        super().__init__()\n",
    "        self.n_layer = n_layer\n",
    "        self.jk = jk  # jumping-knowledge\n",
    "        self.residual = residual  # residual/skip connections\n",
    "        self.atom_encoder = AtomEncoder(emb_dim=emb_dim)  # encode node features\n",
    "        self.convs = nn.ModuleList([GINConv(emb_dim) for _ in range(n_layer)])\n",
    "        norm = nn.BatchNorm1d if norm == 'batch' else LayerNorm\n",
    "        self.bns = nn.ModuleList([norm(emb_dim) for _ in range(n_layer)])\n",
    "        self.pool = global_add_pool\n",
    "\n",
    "    def forward(self, data: Batch) -> torch.Tensor:\n",
    "        x, edge_index, edge_attr, batch = data.x, data.edge_index, data.edge_attr, data.batch\n",
    "        # node embeddings\n",
    "        hs = [self.atom_encoder(x)]\n",
    "        for layer in range(self.n_layer):\n",
    "            h = self.convs[layer](hs[layer], edge_index, edge_attr)\n",
    "            h = self.bns[layer](h)\n",
    "            if layer < self.n_layer - 1:\n",
    "                h = F.relu(h)\n",
    "            if self.residual:\n",
    "                h += hs[layer]\n",
    "            hs += [h]\n",
    "        # graph embedding and prediction\n",
    "        if self.jk:\n",
    "            h = torch.cat([h for h in hs], -1)\n",
    "        h_out = self.pool(h, batch)\n",
    "        return h_out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "84a30106",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:alibi_detect.cd.utils:Input shape could not be inferred. If alibi_detect.models.tensorflow.embedding.TransformerEmbedding is used as preprocessing step, a saved detector cannot be reinitialized.\n"
     ]
    }
   ],
   "source": [
    "from alibi_detect.cd import MMDDrift\n",
    "\n",
    "enc = Encoder(n_layer=1).to(device)\n",
    "preprocess_fn = partial(predict_batch, model=enc, device=device, preprocess_fn=batch_fn, batch_size=32)\n",
    "dd = MMDDrift(x_ref, backend='pytorch', p_val=.05, n_permutations=1000, preprocess_fn=preprocess_fn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "28169267",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "H0 ROC-AUC: 0.763 -- loss: 0.208\n",
      "Drift? No!\n",
      "p-value: 0.925\n",
      "\n",
      "test sample ROC-AUC: 0.736 -- loss: 0.150\n",
      "Drift? Yes!\n",
      "p-value: 0.000\n",
      "\n",
      "imbalanced sample ROC-AUC: 0.749 -- loss: 0.940\n",
      "Drift? Yes!\n",
      "p-value: 0.000\n",
      "\n"
     ]
    }
   ],
   "source": [
    "make_predictions(dd, xs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0af91e59",
   "metadata": {},
   "source": [
    "### Input data drift using a learned kernel\n",
    "\n",
    "Instead of applying the MMD detector on the pooling output of a randomly initialized GNN encoder, we use the [Learned Kernel detector](https://docs.seldon.io/projects/alibi-detect/en/stable/cd/methods/learnedkerneldrift.html) which trains the encoder and kernel on part of the data to maximise an estimate of the detector's test power. Once the kernel is learned a permutation test is performed in the usual way on the value of the MMD on the held out test set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "e7aea114",
   "metadata": {},
   "outputs": [],
   "source": [
    "from alibi_detect.cd import LearnedKernelDrift\n",
    "from alibi_detect.utils.pytorch import DeepKernel\n",
    "\n",
    "kernel = DeepKernel(enc, kernel_b=None)  # use the already defined random encoder in the deep kernel\n",
    "dd = LearnedKernelDrift(x_ref, kernel, backend='pytorch', p_val=.05, dataloader=DataLoader, \n",
    "                        preprocess_batch_fn=batch_fn, epochs=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "96f9cc64",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "H0 ROC-AUC: 0.763 -- loss: 0.208\n",
      "Drift? No!\n",
      "p-value: 0.630\n",
      "\n",
      "test sample ROC-AUC: 0.736 -- loss: 0.150\n",
      "Drift? Yes!\n",
      "p-value: 0.000\n",
      "\n",
      "imbalanced sample ROC-AUC: 0.749 -- loss: 0.940\n",
      "Drift? Yes!\n",
      "p-value: 0.000\n",
      "\n"
     ]
    }
   ],
   "source": [
    "make_predictions(dd, xs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d73d6be",
   "metadata": {},
   "source": [
    "Since the molecular scaffolds are different across the train, validation and test sets, we expect that this type of data shift is picked up in the input data (technically not the input but the graph embedding).\n",
    "\n",
    "### Drift on graph statistics\n",
    "\n",
    "We could also compute graph-level statistics such as the number of nodes, edges and clustering coefficient and detect drift on those statistics using the Kolmogorov-Smirnov test with multivariate correction (e.g. [Bonferroni](https://en.wikipedia.org/wiki/Bonferroni_correction)). First we define a preprocessing step to extract the summary statistics from the graphs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "a4d08875",
   "metadata": {},
   "outputs": [],
   "source": [
    "# return number of nodes, edges and average clustering coefficient per graph\n",
    "def graph_stats(data: List[Data]) -> np.ndarray:\n",
    "    num_nodes = np.array([d.num_nodes for d in data])\n",
    "    num_edges = np.array([d.num_edges for d in data])\n",
    "    c = np.array([np.array(list(clustering(to_networkx(d)).values())).mean() for d in data])\n",
    "    return np.concatenate([num_nodes[:, None], num_edges[:, None], c[:, None]], axis=-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "7ee6d196",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:alibi_detect.cd.utils:Input shape could not be inferred. If alibi_detect.models.tensorflow.embedding.TransformerEmbedding is used as preprocessing step, a saved detector cannot be reinitialized.\n"
     ]
    }
   ],
   "source": [
    "dd = KSDrift(x_ref, p_val=.05, preprocess_fn=graph_stats)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "ee7ceea2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "H0 ROC-AUC: 0.763 -- loss: 0.208\n",
      "Drift? No!\n",
      "p-value: [0.9884282  0.79581064 1.        ]\n",
      "\n",
      "test sample ROC-AUC: 0.736 -- loss: 0.150\n",
      "Drift? No!\n",
      "p-value: [0.94992185 0.23822938 0.52494323]\n",
      "\n",
      "imbalanced sample ROC-AUC: 0.749 -- loss: 0.940\n",
      "Drift? Yes!\n",
      "p-value: [5.1066105e-05 1.1882804e-05 9.5120013e-01]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "make_predictions(dd, xs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd867e31",
   "metadata": {},
   "source": [
    "The 3 returned p-values correspond to respectively the p-values for the number of nodes, edges and clustering coefficient. We already saw in the EDA that the distributions of the node, edge and clustering coefficients look similar across the train, validation and test sets except for the tails. This is confirmed by running the drift detector on the graph statistics which cannot seem to pick up on the differences in molecular scaffolds between the datasets, unless we heavily oversample from the minority class where the number of nodes and edges but not the clustering coefficient significantly differ."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
